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1.  INTRODUCTION 


1.1  This  Work  vs.  ELatsambalos  (1981) 

This  work  is  concerned  with  the  computation  of  the  gravity 
disturbance  vector  external  to  the  earth,  from  given  gravity  anomaly 
data  referred  to  the  earth's  surface.  The  central  problem  is  how  to 
rigorously  account  for  the  irregular  shape  of  the  topographic  surface  to 
which  the  data  refer.  This  work  is  a  continuation  of  the  work  of 
Katsambalos  (1981)  entitled  "Simulation  Studies  on  the  Computation  of 
the  Gravity  Vector  in  Space  from  Surface  Data  Considering  the 
Topography  of  the  Earth". 

There  are  a  number  of  differences  between  the  present  work  and 
that  of  Katsambalos.  One  difference  is  the  emphasis  placed  in  the 
present  study  upon  the  use  of  suitable  models  for  the  separate  modeling 
of  different  frequency  ranges  of  the  total  disturbance  vector  signal.  To 
this  end,  greater  attention  than  that  in  Katsambalos’  Section  9.3  is  now 
given  to  the  use  of  the  spherical  harmonic  representation  to  model  the 
low  frequency  components  of  the  signal.  And,  as  an  entirely  new 
addition  to  the  work  of  Katsambalos,  a  study  is  now'  made  of  modeling 
the  very  high  frequency  components  of  the  signal  as  the  integrated 
gravitational  effects  of  certain  shallow  topographic  masses  of  assumed 
constant  density.  The  representation  by  topographic  mass  effects,  in 
which  the  input  consists  of  detailed  terrain  elevation  data,  is  advocated 
in  the  works  of  Tscherning  and  Porsberg  (see  Tscherning  (1979), 
Forsberg  and  Tscherning  (1981),  Tscherning  and  Porsberg  (1983), 
Porsberg  (1984)).  The  residual  signal  field,  left  after  subtracting  from 
the  total  field  the  fields  already  represented  by  spherical  harmonics  and 
topographic  mass  effects,  is  a  relatively  smooth  and  low  energy  field. 

As  suitable  models  for  representing  the  residual  field,  a  study  is 
made  here  of  the  classical  integral  model,  the  so~called  Dirac  approach, 
and  the  least  squares  collocation  approach,  which  are  models  studied 
also  in  Katsambalos  (ibid).  In  contrast  to  Katsambalos,  the  present 
study  relies  upon  the  use  of  the  classical  integral  model,  but  with  mean 
topography  accounted  for,  as  one  of  viable  complementary  models  that 
can  be  operationally  used.  For  the  Dirac  and  least  squares  collocation 
approaches  the  present  study  gives  an  expanded  presentation  over  that 
of  Katsambalos,  in  order  to  be  able  to  numerically  experiment  on 
different  possible  models  that  actually  fall  under  the  heading  of  being 
Dirac  or  least  squares  collocation  approaches.  Specifically,  point  mass 
modeling  and  point  dipole  modeling  are  experimented  on  under  the 
heading  of  being  Dirac  approaches,  in  addition  to  experimenting  on  the 


1 


V.v 


.■.V. 

•* 


.•V 


» •• 


2 


usual  modeling  by  gravity  anomaly  impulses  on  an  internal  sphere 
(Kataambalos,  ibid,  Chapter  6).  Under  the  least  squares  collocation 
approach,  experiments  are  performed  using  empirical  covariance 
functions  that  are  based  on  white  noise  gravity  anomaly  and  white  noise 
disturbing  potential  distributions  on  an  internal  sphere.  Finally,  the 
present  study  has  taken  into  account  the  questions  raised  by 
Bjerhammar  and  Sjoberg  (1982,  private  communications)  and  by 
Tscherning  (1983a)  about  Katsambalos*  conclusions  on  the  Dirac  and  least 
squares  collocation  approaches. 

For  numerical  experimentations  the  present  study  has  now  used  real, 
as  opposed  to  simulated,  gravity  and  terrain  elevation  data.  This  allows 
for  the  verification  and  evaluation  of  various  theoretical  models  as 
applied  to  the  actual  variations  present  in  the  earth’s  gravity  field.  A 
detailed  study  of  the  spectral  characteristics  of  the  spatial  disturbance 
vector  signal,  as  well  as  its  response  to  gravity  anomaly  data  of  varying 
resolution  and  distance  away  from  the  computation  point,  is  given 
(Chapter  3)  to  help  in  the  design  of  models  and  numerical  experiments. 

It  should  be  mentioned  that  the  Green's  approach  studied  also  in 
Katsambalos  is  not  at  ail  used  in  the  present  study.  The  reason  is  that 
the  said  approach  requires  as  data  gravity  disturbances  and  disturbing 
potentials  on  the  earth’s  surface,  the  elevations  of  the  surface  points, 
and  the  North-South  and  East-West  components  of  the  surface 
inclinations.  These  data  are  very  difficult  or  impossible  to  obtain  with 
sufficient  accuracy  and  density  in  an  operational  environment. 


1.2  Preliminaries 

We  use  the  following  conventional  notations  of  gravimetric  geodesy: 

^  gravity  vector  of  the  earth,  referred  to  simply  as  gravity  vector 

^  gravity  vector  of  the  reference  ellipsoid,  referred  to  as  normal 

gravity  vector 

g  magnitude  of  referred  to  as  gravity 

y  magnitude  of  referred  to  as  normal  gravity 

W  potential  of  t  (i«e.,  t  -  grad  W),  referred  to  as  gravity  potential 

U  potential  of  ^  (i.e.,  ^  =  grad  U),  referred  to  as  normal  gravity 
potential. 

In  the  above  we  are  working  in  a  body-fixed  coordinate  system  (i.e.,  one 
rotating  with  the  earth  or  ellipsoid),  and  the  gravity  vector  is  the 
resultant  of  the  body’s  gravitational  force  and  the  centrifugal  force  of 
the  system’s  rotation. 
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Fundamental  to  most  gravimetric  analyses  is  the  disturbing  potential 
Ti  defined  as 


T  =  W  -  U. 


(1.1) 


The  force  associated  with  T  is  the  gravity  disturbance  vector  'i,  i.e., 


^  =  grad  T. 


(1.2) 


Based  on  the  above  definitions,  ^  can  also  be  written  in  the  familiar 
form: 


(1.3) 


Equation  (1.3)  forms  a  basis  for  modeling  the  gravity  vector  as 
the  sum  of  ^  and  t.  The  sum  is  dominated  by  the  normal  gravity  vector 
which  is  rigorously  computable  (see,  e.g.,  Heiskanen  and  Moritz,  1967, 
sec.  6-2).  For  the  spatial  gravity  disturbance  vector  no  closed 
computational  formulas  exist,  but  its  accurate  modeling  from 
observational  data  on  the  earth’s  surface  is  precisely  the  subject  of  the 
present  study. 

Basically,  the  modeling  of  "S  proceeds  through  its  relation  (1.2)  to 
the  fundamental  function  T.  In  turn  the  function  T  can  be  linked  to 
quantities  obtainable  from  observations,  most  commonly  the  gravity 
anomaly.  Indeed,  in  this  report  we  use  the  gravity  anomaly  as 
fundamental  data  for  our  modeling  procedures.  In  the  following  we 
detail  the  precise  meaning  of  the  quantity  referred  to  as  gravity 
anomaly,  and  discuss  the  relationship  between  gravity  anomaly  and 
disturbing  potential  T. 

We  have  the  following  conventional  definition  of  gravity  anomaly  at  a 
point  P  (Heiskanen  and  Moritz,  1967,  pp.  83,  91;  Moritz,  1980,  p.  353; 
Jekeli,  1981,  pp.  39,  120;  Moritz,  1983,  pp.  4-7): 


4gp  =  ip  -  ro- 


(1.4) 


That  is,  the  gravity  anomaly  at  P  is  the  difference:  gravity  at  P  minus 
normal  gravity  at  some  point  Q.  Point  Q  is  called  the  normal  point  of  P, 
and  is  established  such  that  (a)  the  normal  gravity  potential  U  at  Q  is 
equal  to  the  actual  gravity  potential  W  at  P  (see  Figure  1): 


U(Q)  =  W(P), 


. 


(1.5) 


Normal  Plumb 
Line 


U  *  U(Q)  =  W(P) 
(spheropotential  surface) 


Figure  1.  Construction  of  the  Gravity  Anomaly. 


and  that  (b)  P  and  Q  lie  on  the  same  plumb  line  of  the  normal  gravity 
field.  The  gravity  anomaly  (1.4)  is  distinct  from  the  gravity  anomaly 
vector,  defined  as  (Heiskanen  and  Moritz,  1967,  p.  83): 


AJp  =  ip  ~  fg. 


(1.6) 


Jekeli  (1981,  Chapter  4)  shows  that  the  definition  (1.4)  leads  to  the 
following  rigorous  relation  between  the  gravity  anomaly  and  the 
disturbing  potential  at  a  point: 


Agp  =  - 


1  *Tp 
7p  'hp 


Tp 


+  «p  +  0(e*), 


(1.7) 


with 


*p  = 


(1.8) 


and  where 


•  derivative  along  the  plumb  line  of  the  actual  gravity  field 

*  derivative  along  the  plumb  line  of  the  normal  gravity  field; 

4h  this  is  not  distinguished  from  the  derivative  along  the 

straight  ellipsoidal  normal  (see  ibid.,  p.ll9) 

0(e*)  indicates  that  the  omitted  terms  are  on  the  order  of  e*  Ag<0.01 
mgal  (e*  =  square  of  eccentricity  *  0.0067;  and  dg  has  a 
maximal  value  of  300  mgal). 

The  term  ep  of  (1.8)  is  dependent  on  the  deflections  of  the  vertical 
(ibid.,  eqs.  (4.8),  (4.22)),  and  has  a  magnitude  of  less  than  0.34  mgal. 
Retaining  only  the  first  two  terms  and  dropping  the  subscripts  in  (1.7), 
we  obtain: 


ill 

y  *h 


T. 


(1.9) 


Equation  (1.9)  is  exactly  the  same  as  the  boundary  condition  that 
results  from  conventional  linearizations  of  Molodensky’s  boundary  value 
problem  (see  e.g.,  Molodenskii  et  al.,  1962,  chap.  V,  or  Heiskanen  and 
Moritz,  1967,  sec.  8-5).  Moritz  (1980,  p.336)  writes  that  such 

lineeu'izations  are  "practically  sufficiently  accurate  but  not  completely 
rigorous".  In  the  rigorous  linearization  of  Krarup  (1973)  (see  Moritz, 
ibid.,  pp.  337-349)  the  resulting  boundary  condition  is  of  the  form  (1.9) 
also,  but  with  Ag  interpreted  as  the  component  of  the  gravity  anomaly 
vector  (see  (1.6))  in  the  downward  direction  of  the  isozenithal  and  with 
the  derivatives  taken  along  the  direction  of  the  isozenithal.  For 
practical  modeling  problems,  which  is  our  interest  in  the  present  report, 
it  is  convenient  to  use  (1.9)  in  our  original  meaning  that  4g  is  the 
gravity  anomaly  (1.4)  and  that  derivatives  are  taken  along  the  ellipsoidal 
normal.  Moritz  (ibid.,  pp.  352-353)  has  termed  (1.9)  as  the  "practical 
boundary  condition". 

Under  the  so-called  spherical  approximation  (Heiskanen  and  Moritz, 
1967,  pp.  87-88),  the  normal  derivative  is  approximated  by  the  radial 
derivative. 


9 

»h 


9 

Tp’ 


(1.10) 


and  the  normal  gravity  is  approximated  by  the  attraction  of  a 
homogeneous  sphere, 


with  kM  being  the  geocentric  gravitational  constant.  Under  (1.10)  and 
(1.11),  (1.9)  transforms  to  the  usual  spherical  approximation: 


Ag  =  - 


II 

*r 


(1.12) 


From  Jekeli  (ibid.,  pp.  122-123)  a  more  rigorous  transformation  of  (1.9) 
neglecting  only  terms  of  0(e*)  can  be  written  as  follows: 


aT  2 

Ag=-— --T-e»  sin* 
+  [sjj  fr  P,(sin*)  - 


To  review  standard  notations, 

r,  ♦,  X  geocentric  radius,  geocentric  latitude,  and  longitude  of  the  point 
to  which  the  Ag  and  T  refer 

e  eccentricity  of  the  reference  ellipsoid 

Jj  second  degree  zonal  harmonic  of  the  normal  gravity  field 

P2  second  degree  Legendre  polynomial 

a  semi-major  axis  of  reference  ellipsoid 

kM  geocentric  gravitational  constant 

rotational  speed  of  the  earth. 

Equation  (1.13)  is  of  the  form  of  (1.12),  but  with  two  terms  on  the  order 
of  e*  added  to  correct  for  the  effect  of  spherical  approximation  in 
(1.12).  To  see  that  the  last  term  of  (1.13)  is  O(e’),  refer  to  Equations 
(2.14)  and  (2.15)  of  Chapter  2.  The  first  correction  term  of  (1.13) 
corrects  the  first  term  of  (1.12),  and  the  second  correction  term 
corrects  the  second  term  of  (1.12).  If  terms  on  the  order  of  the  earth’s 
flattening  f  (f  *  e*/2  »  0.003)  are  neglected,  then  (1.13)  reduces  to 
(1.12). 

Indeed,  the  formal  meaning  of  spherical  approximation  is  the 
neglection  of  terms  on  the  order  of  f  in  equations  such  as  (1.13) 
relating  quantities  of  the  anomalous  gravity  field.  The  anomalous 


-  ^T 

cos*  - :•  + 

ra^ 


3  ^  (1  -  sin»*)]T. 


(1.13) 


we  have: 


gravity  field  is  the  one  associated  with  the  disturbing  potential  T.  A 
useful  geometric  meaning  of  spherical  approximation  is  given  in  Moritz 
(1980)  pp.  351-352).  This  meaning  is  in  terms  of  a  spatial  mapping  of 
pointSi  as  shown  in  Figure  2. 

It  is  shown  that  a  point  P  with  geocentric  coordinates  (r,  X)  is 
mapped  onto  a  point  P'  with  geocentric  coordinates  (R  *  h,  X),  where 
R  is  a  mean  earth  radius  (usually.  R  =  6371  km),  h  is  the  height  of 
point  P  above  the  reference  ellipsoid,  and  4  is  the  geodetic  latitude  of 
P.  The  spherical  approximation  consists  in  using  P'  instead  of  P  in  all 
calculations.  Formally,  this  means  using  (R  *  h,  ^)  as  geocentric 
(radius,  latitude)  in  all  calculations. 

Finally,  below  we  outline  the  procedure  by  which  the  gravity 
anomaly  as  defined  in  (1.4)  can  be  obtained  from  observational  data  on 
the  earth’s  surface.  More  details  can  be  found  in  Rapp  (1984,  pp.  3-7). 
The  procedures  result  in  gravity  anomalies  that  refer  to  the  earth’s 
surface,  and  these  anomalies  are  the  fundamental  data  for  our  modeling 
procedures  in  this  report.  We  have  the  following  steps: 


P'(R+h,  4>,  x) 


ELLIPSOID  SPHERE 


Figure  2.  Mapping  Under  the  Spherical  Approximation 


a)  gravity  gp  is  measured  on  the  earth's  surface 

b)  by  leveling  procedures  the  potential  difference  denoted  by  Cpj  is 
obtainedt  where 


Cpi  =  Wp  -  Woi  (1.14) 

with 

Wp  gravity  potential  at  P 

Woi  gravity  potential  on  the  iiil  equipotential  surface  of  the 
earth’s  gravity  field  being  used  as  reference  surface  of 
leveling.  The  reference  surface  of  leveling  is  also  referred 
to  as  "height  datum"  or  "vertical  datum".  The  index  i  is 
used  to  account  for  the  realistic  case  of  several  different 
height  datums  being  used  in  different  parts  of  the  world, 
i.e.,  there  is  no  unique  world  vertical  datum  at  present. 

c)  A  "normal  height"  Hj*  (Heiskanen  and  Moritz,  1967,  sec.  4-5)  can  be 
obtained  from: 


(1.15) 


where  y  is  an  average  value  of  normal  gravity  along  the  plumb  line, 

d)  A  normal  gravity  is  computed: 


ri  =  r(Hi*,  tp,  Xp),  (1.16) 

i.e.,  as  the  normal  gravity  at  geodetic  coordinates  (Hi»,  4p,  Xp).  In 
practice  ,  it  is  satisfactory  here  to  use  instead  of  Hi*  the 
orthometric  height  (Cruz  and  Laskowski,  1984,  eq.  (3.1.4)). 

e)  A  gravity  anomaly  is  computed: 


4gi  =  gp  -  ri- 


(1.17) 


f)  The  Agj  of  (1.17)  should  then  be  corrected  to  yield  the  gravity 
anomaly  Ag  defined  by  (1.4).  We  first  define  (Rapp,  1984,  eqs.  (13), 


4W  =  Wo  -  Uo 


AWoi  =  Wo  Wo if 


(1.18) 


(1.19) 


where 

Wo  gravity  potential  on  an  assumed  world  vertical  datum  (e.g.,  the 
geoid ) 

Uo  normal  gravity  potential  on  the  surface  of  the  equipotential 
reference  ellipsoid. 

Theni  the  Ag  of  (1.4)  can  be  obtained  from  the  Agi  of  (1.17)  as 
follows  (see  Rapp,  ibid.,  eq.  (29)): 


Ag  =  Agi  +  - 


-  AW. 
r 


(1.20) 


g)  Theoretical  developments  in  gravimetric  geodesy  have  been 
considerably  simplified  by  assuming  the  disburbing  potential  T  to  be 
a  harmonic  function  outside  the  earth’s  attracting  masses.  To 
satisfy  this  assumption  the  mass  of  the  atmosphere  must  be 
computationally  removed  from  the  problem.  As  it  is,  for  example,  the 
force  of  "gravity"  measured  in  step  (a)  includes  the  attraction  of 
the  atmosphere  at  the  point  P.  The  removal  of  the  atmosphere  from 
the  main  modeling  problem,  and  its  subsequent  restoration  in  the 
final  result,  can  be  done  in  a  simple  way  because  of  the  smallness 
of  atmospheric  effects.  The  theory  is  discussed  in  Moritz  (1980,  pp. 
422-425).  As  a  practical  procedure  for  removing  atmospheric  effects 
from  gravity  anomaly  data,  the  mass  of  the  atmosphere  is  first 
added  to  the  mass  of  the  earth  to  form  the  defined  mass  of  the 
reference  ellipsoid.  This  affects  the  value  of  in  (1.17).  Then  an 
"atmospheric  correction"  6g/^  is  added  to  the  computed  gravity 
anomaly  Agi  of  (1.17)  to  get  the  corrected  anomaly  Agj*: 


Agi*  =  Agi  +  <5gA. 


(1.21) 


Values  of  (<1  mgal)  are  tabulated  as  a  function  of  height  h  in 
lAG  (1970,  pp.  72-73).  The  Agi*  instead  of  igi  is  the  value  to  use  in 
(1.20)  to  yield  the  Ag  as  defined  by  (1.4),  under  the  assumption 
that  there  is  no  atmosphere. 

h)  Another  simplifying  procedure  in  gravimetric  computations  is  the 
removal  of  the  zero  and  first  degree  spherical  harmonic  terms  of  the 
disturbing  potential  T  from  the  main  modeling  problem.  These  terms 
can  be  modeled  using  a  separate  set  of  procedures  (see,  e.g., 
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Heiskanen  and  Moritz,  1967,  pp.  98-99,  107).  The  first  decree  terms 
do  not  specifically  affect  the  gravity  anomaly  Ag,  which  is  related  to 
T  through  (1.12).  The  effect  of  the  zero  degree  term  can  be 
removed  from  gravity  anomaly  data  Agi,  to  get  the  corrected  value 
Agi',  as  follows  (see  Rapp,  1984,  p.8): 


=  Agi  +  ^  , 


(1.22) 


where  AkM  =  kM  -  kMg;  kMg  is  the  gravitational  constant  of  the 
reference  ellipsoid. 

i)  Summarizing,  the  gravity  anomaly  that  satisfies  (1.4)  under  the 
simplification  of  steps  (g)  and  (h)  is: 


(1.23) 


Ag  =  Agi  +  Ag*  +  “AWoi  -  ^AW  + 


where  the  Agi  is  from  (1.17). 


1.3  Scope  of  This  Study 

Gravity  anomalies  that  refer  to  the  earth’s  surface  will  be  the 
fundamental  data  for  our  modeling  procedures  in  this  report.  The 
topography  itself  is  defined  through  the  ellipsoidal  heights  of  surface 
points.  Such  heights  are,  for  the  purposes  of  analyzing  the  anomalous 
gravity  field,  not  distinguished  from  available  orthometric  heights  or 
some  approximate  heights  from  topographic  maps  or  digital  terrain 
models  (Moritz,  1980,  p.  357). 

First,  we  account  for  the  fact  that  spherical  harmonic  series 

representation  has  turned  out  to  be  a  very  efficient  representation  of 
the  anomalous  gravity  field  up  to  some  limited  resolution  (see,  e.g., 

Tscherning,  1983b).  Therefore,  in  Chapter  2  we  discuss  the 
re-parameterization  of  the  long-wavelength  part  of  the  anomalous 

gravity  field  from  the  given  surface  gravity  anomalies  to  a  set  of 
spherical  harmonic  coefficients  of  the  earth’s  disturbing  potential. 

Currently  envisioned  maximum  degree  of  spherical  harmonic  expansions 
is  around  n  =  360,  corresponding  roughly  to  a  30'x30'  resolution. 
Theoretical  and  practical  problems  associated  with  the  determination  of 
the  coefficients  in  such  high  degree  expansions,  using  a  combination  of 
coefficients  derived  from  satellite  motion  analysis  and  of  coefficients 
implied  by  terrestrial  gravity  anomaly  information,  are  recently  reviewed 
and  discussed  in  Rapp  (1984).  Our  specific  concern  in  Chapter  2  is  the 
computation  of  spherical  harmonic  coefficients  from  surface  gravity 
anomaly  data,  in  a  manner  that  takes  into  account  the  problems  of 
ellipticity  and  topography  of  the  earth.  These  problems  are  also 
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discussed  in  Rapp  (ibid.),  with  references  to  the  treatment  of  ellipticity 
by  Pellinen  (1982).  Our  contribution  via  Chapter  2  is  in  providing  a 
geometric  interpretation  to  the  so-called  ellipsoidal  correction  terms 
derived  in  Pellinen  (ibid.),  by  using  elementary  Taylor  series 
considerations  in  re-deriving  those  terms.  Using  this  geometric 
interpretation  in  combination  with  the  treatment  of  topography  in  Rapp 
(1984)  we  form  recommendations  on  the  computation  of  harmonic 
coefficients  from  surface  gravity  anomaly  data.  These  recommendations 
have  some  differences  with  those  of  Rapp  (ibid.). 

The  spherical  harmonic  representation  is  feasible  only  up  to  some 
limited  resolution,  currently  on  the  order  corresponding  to  a  maximum 
harmonic  degree  of  about  n  =  360.  Problems  that  increase  in  importance 
with  unlimited  increase  in  the  degree  of  spherical  harmonic 
representation  include  (a)  the  effects  of  series  divergence  below  the 
smallest  sphere  bounding  the  earth  (see,  e.g.,  Jekeli,  1981;  Colombo, 
1982);  (b)  the  computer  storage  and  time  requirements  for  handling  a 
very  large  number  of  coefficients  (see  Tscherning  et  al.,  1983);  and  (c) 
the  determination  of  the  coefficients  of  the  series,  causing  problems  on 
computer  requirements,  numerical  quadratures  and  the  implementation 
of  a  relation  between  coefficients  and  gravity  anomaly  data  that 
accounts  for  the  earth's  ellipticity  and  topography  (see  Rapp,  1984). 

Therefore,  to  complete  the  modeling  of  the  gravity  field  up  to  very 

detailed  resolutions,  the  globally  valid  spherical  harmonic  representation 
must  be  complemented  by  locally  valid  solutions. 

Before  embarking  on  a  detailed  study  of  various  locally  valid  models, 
in  Chapter  3  we  first  present  a  familiarization  study  of  our  signal  of 
interest,  namely,  the  gravity  disturbance  vector  in  space.  The 

information  in  Chapter  3  will  help  in  the  design  of  solutions  and 

numerical  experiments.  Spherical  earth  formulas  for  spherical  harmonic 
and  integral  representations  of  the  disturbance  vector  are  first 

reviewed  in  Section  3.1.  The  neglect  of  the  earth  ellipticity  and 

topography  should  not  affect  the  conclusions  of  the  chapter,  and  at  the 
same  time  allows  for  the  use  of  spectral  analysis.  In  Section  3.2, 

spectral  analysis  yields  the  power  of  the  signal  within  a  given 
frequency  range,  as  a  function  of  altitude  in  space.  Knowledge  of  such 
power  distributions  is  useful  for  determining  the  model  resolutions 
required  for  representing  all  significant  signal  energies  at  a  given 
altitude.  In  Section  3.3  the  so-called  truncation  theory  yields  the 

sensitivity  of  the  signal  at  a  given  altitude,  to  gravity  anomaly  data  of 
certain  resolution  and  distance  away  from  the  computation  point.  Such 
information  is  useful  for  minimizing  the  required  data  cap  size  for  a 
given  resolution  of  gravity  anomaly  data.  In  Section  3.3  we  also  examine 
the  use  of  different  truncation  theories,  namely,  the  unmodified 
Molodensky  truncation,  Meissl  truncation,  Wong/Gore  truncation,  and 
Improved  Molodensky  truncation  as  applied  to  disturbance  vector 
computations  at  altitude. 

In  Chapter  4  we  begin  to  look  at  the  modeling  of  the  entire 
frequency  range  of  the  disturbance  vector  signal,  from  low  to  high 
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frequencies.  In  Section  4.1  some  general  considerations  are  first  given 
to  help  understand  the  effect  of  the  non-spherical  shape  of  the 
boundary  surface  (i.e.,  the  topography  to  which  the  given  gravity 
anomaly  data  refer)  on  the  modeling  of  the  external  disturbance  vector 
field.  Then  in  Section  4.2  a  composite  model  is  presented  that  in 
principle  models  the  total  disturbing  potential  T  in  space  as  the  sum  of 
three  components: 

T  =  TS  +  TO  +  T*.  (1.24) 


where 


T®  long  wavelength  component,  generated  from  the  spherical 
harmonic  representation; 

short  wavelength  component,  generated  as  the  potential  of 
certain  shallow  topographic  masses  of  assumed  constant  density; 

TO  medium  wavelength  residual  component,  equal  to  whatever  is  the 
difference  (T  -  T®  -  T*).  The  TO  is  modeled  in  Section  4.2.2 
using  the  classical  direct  integration  of  gravity  anomaly  data, 
with  mean  topography  accounted  for. 

The  gravity  anomalies  AgO  that  are  integrated  to  yield  tO  are  expressed 
through  (1.12): 


Ago  =  -  ^  ^  TO.  (1.25) 


which  by  definition  of  TO  becomes: 


or 

AgO  =  Ag  -  Ag®  -  Ag‘, 


(1.26) 


(1.27) 


i.e.,  the  AgO  is  the  residual  anomaly  left  after  subtracting  from  the 
original  anomaly  Ag  the  anomalies  Ag®  implied  by  T®  and  Ag^  implied  by 
T^.  Note  that  AgO  continues  to  refer  to  the  earth’s  surface. 


The  components  of  the  disturbance  vector  corresponding  to  T®,  T’^i 
and  are: 


(Ifs,  ^0,  t*)  =  gradCTS,  to,  T*),  (1.28) 


and  the  total  disturbance  vector  is  then  modeled  as: 


^  +  ^0  +  2! 


In  Section  4.3  numerical  experiments  are  conducted  on  the  complete 
modeling  of  "i  using  (1.29).  In  all  experiments  real  gravity  anomaly  and 
elevation  data  are  used  over  a  7*x9*  mountainous  test  area  in  New 

Mexico.  Of  special  interest  are  the  experiments  that  give  a  feeling  for 

how  well  the  l^o  ^an  absorb  the  resulting  unmodeled  part  caused  by  the 
omission  of  either  it®  or  in  (1.29).  The  ability  of  ijo  absorb  7^®  is 

mainly  related  to  the  size  of  the  anomaly  data  cap  used  in  The 

ability  of  to  absorb  is  related  to  the  modeling  error  introduced  by 
accounting  only  for  a  mean  topography,  and  not  for  the  full  variations 
of  the  topograpy,  in  computing  from  the  surface  data  Ag“. 

The  solution  expressed  by  (1.24)  may  be  called  a  continuous 
approach,  since  the  underlying  theories  assume  a  continuous  coverage  of 
the  earth’s  surface  by  gravity  anomaly  data.  Under  this  continuous 
approach  the  integration  of  anomaly  data  for  the  computation  of  the 
residual  potential  T°  must  strictly  be  performed  along  the  earth’s 
surface  where  the  data  refer.  Molodensky  and  later  Brovar  have 
written  the  solution  in  terms  of  a  series  of  integrals  (see  Moritz,  1980, 
secs.  43  and  44,  with  reference  to  Molodenskii,  et  al.,  1962  and  Brovar, 
1964).  However,  there  is  an  equivalent  solution  which  is  convenient  for 
conceptualizations  and  discussions,  because  it  is  full  of  physical 
meaning.  Details  of  this  solution  are  given  in  Heiskanen  and  Moritz 
(1967,  sec.  8-10  and  pp.  324-325,  with  reference  to  Bjerhammar,  1964) 
and  in  Moritz  (1980,  secs.  45,  46  and  pp.  419-420).  The  solution  consists 
of  (a)  analytical  continuation  of  surface  anomaly  data  to  a  level  surface, 
followed  by  (b)  classical  integration  along  the  level  surface  as 
approximated  by  a  sphere  (see  Figure  2). 

In  practice,  the  computation  of  analytically  continued  anomalies 
poses  rather  severe  data  density  and  accuracy  requirements  especially 
in  areas  of  rough  gravity  field  (see,  e.g.,  Noe,  1980).  Therefore,  the 
motivation  of  the  form  (1.24)  is  to  make  the  field  T*^  as  smooth  as 
possible,  by  "moving"  the  high  frequency  components  of  T  away  from  T° 
into  T^.  With  T'^  smooth,  the  difference  between  the  surface  sinomalies 
Ag*^  of  (1.27)  and  the  corresponding  analytically  continued  anomalies  can 
be  so  small  as  to  be  negligible.  This  way  we  avoid  the  specific 
computation  of  analytically  continued  anomalies  by  simply  equating  these 
anomalies  to  the  surface  anomalies.  Numerical  investigations  related  to 
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these  ideas  are  also  included  in  Chapter  4. 

A  more  rigorous  way  to  avoid  the  above  analytical  continuation 
computations  of  the  continuous  approach  would  be  to  perform  collocation 
in  space  (see  Moritz,  1983,  p.  31).  The  theory  is  reviewed  and  numerical 
investigations  performed  in  Chapter  5  of  this  study.  Two  approaches  to 
collocation  are  studied  here,  namely,  the  so-called  Dirac  approach  and 
the  least  squares  collocation  approach.  Collocation  is  a  discrete 
approach,  i.e.,  the  theory  assumes  the  data  to  be  given  only  at  a  finite 
number  of  discrete  points  on  the  earth's  surface.  The  collocation 
solution  yields  an  approximation  to  the  disturbing  potential,  that  (a)  is 
harmonic  down  to  an  internal  sphere  (spherical  approximation)  completely 
embedded  within  the  earth  and  (b)  satisfies  the  data  at  the  given 
points.  For  feasible  economy  and  convergence  of  solution,  it  is  still 
desirable  to  use  the  complementary  models  T^  and  T^  as  much  as 
possible,  as  in  (1.24).  In  this  case  we  have  the  model: 


T  =  T®  +  TC  +  (1.30) 

where  denotes  the  residual  component 

TC  =  T  -  T*  -  1^,  (1.31) 


as  modeled  by  collocation  procedures  from  the  gravity  anomaly  data 
(see  (1.27)): 


4gC  =  AgO  =  Ag  -  ags  _  Agt.  (1.32) 


Corresponding  to  (1.29)  we  now  have: 


(1.33) 


with 


=  grad  TC.  (1.34) 


As  before,  of  special  interest  in  the  numerical  studies  in  Chapter  5 
would  be  to  see  how  well  can  absorb  the  resulting  unmodeled  part 
caused  by  not  using  "i*-  in  (1.33).  The  discrete  model  would 
expectedly  absorb  better  than  would  the  continuous  model  of 

(1.29),  because  of  the  more  rigorous  way  that  the  topography  is  taken 
into  account  in  the  collocation  procedures.  This  point  is  numerically 
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studied  in  Section  5.4.3. 

Finally,  a  summary  along  with  conclusions  and  recommendations  are 
given  in  Chapter  6. 
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2.  ON  THE  RECOVERY  OF  POTENTIAL  COEFFICIENTS  FROM  SURFACE 

GRAVITY  ANOMALIES 

Various  theoretical  and  operational  aspects  of  solving  the  problem 
have  been  recently  reviewed  by  Rapp  (1984),  discussing  also  the 
combination  of  the  resulting  coefficients  with  those  derived  from  satellite 
motion  analylsis.  In  this  chapter  we  focus  on  one  aspect  of  the 
problem,  namely,  the  development  of  a  mathematical  relation  between 
potential  coefficients  and  surface  gravity  anomaly  data  that  accounts  for 
the  topography  and  ellipticity  of  the  eeurth.  A  feasible  solution  is  to 
first  analytically  continue  the  surface  values  to  values  that  refer  to  a 
sphere,  and  then  use  these  values  on  a  sphere  for  directly  determining 
potential  coefficients  using  the  orthogonality  relationships  of  spherical 
harmonics.  The  question  lies  on  how,  and  to  which  sphere,  the 
analytical  continuation  should  be  done,  and  it  is  to  this  discussion  that 
we  propose  to  contribute  in  this  chapter. 

The  disturbing  potential  T  can  be  represented  in  spherical 
harmonics,  usually  in  the  following  form  which  is  convenient  in 
connection  with  satellite  dynamics  (see  Rapp,  1982b,  (5)): 


T(r,  ♦.  X)  =  CLYnmC^. 


X). 


n=a 


m=-n 


(2.1) 


where 

Cjf„  potential  coefficients.  The  overbar  denotes  full 

normalization,  and  the  asterisk  denotes  the 
subtraction  of  the  potential  coefficients  of  the 
reference  ellipsoid; 


Y„„(4.  X) 


Pniml (C0S4 


cos  mA, 
siniml X, 


m  k  0 

m  <  0 


(2.2) 


Pniml (cos4)  fully  normalized  Legendre  functions  of  degree  n  and 
order  I m I . 


Now  taking  (2.1)  in  conjunction  with  either  (1.12)  or  (1.13),  one  can 
conceive  the  recovery  of  the  coefficients  C5„  from  given  gravity  anomaly 
data.  For  our  purposes  we  will  use  the  more  accurate  (1.13),  especially 
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since  it  turns  out  that  its  more  complicated  form  does  not  add  any 
serious  difficulties  to  the  techniques  of  recovery  that  we  will  discuss. 

First  we  ask,  to  which  sphere  should  the  surface  gravity  anomalies 
be  analytically  continued?  The  answer  is  ultimately  the  equatorial 
sphere  of  radius  r=a,  since  the  coefficients  C*„  as  defined  through  (2.1) 
actually  refer  to  this  sphere.  This  is  the  reason  why,  for  example, 
although  Rapp  (1984,  pp.  17-18)  initially  upward  continued  the  surface 
anomalies  to  the  sphere  of  radius  ~  a  6  km  bounding  the  earth,  the 

factor  (rb/a)*^**  had  to  be  eventually  present  (see  ibid.,  (53))  in  order 

to  downward  continue  the  bounding  sphere  anomalies  back  to  the 

equatorial  sphere  where  the  values  are  directly  needed  for  the 
determination  of  Cj„.  Since  the  sphere  r  =  a  partly  penetrates  the 
earth's  masses,  it  can  be  questioned  whether  the  required  analytical 

continuation  exists  for  the  true  disturbing  potential  T.  However,  here 
we  avoid  this  problem  by  taking  the  view  that  we  are  working  within 
the  context  of  solving  the  Bjerhammar  problem.  That  is,  we  are  using 
the  harmonic  function  T  expressed  by  (2.1)  to  approximate  the  true  T, 
and  such  function  (2.1)  is  indeed  convergent  on  the  sphere  r  =  a  and 
can  be  adjusted  to  agree  with  observational  data  through  the  free 
parameters  Cn„.  The  instability  that  arises  from  the  downward 
continuation  of  surface  data  to  the  embedded  part  of  the  sphere  r  =  a 
is  remedied  by  the  fact  that  we  are  developing  (2.1)  only  up  to  some 
limited  resolution,  expressed  by  truncating  the  series  at  some  maximum 
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Figure  3.  Gravity  Anomalies  on  Different  Boundary  Surfaces. 
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harmonic  degree  Nn,,^  (see  Jekeli,  1981,  pp.  10-11) 

The  next  question  is  how  to  implement  the  analytical  continuation 
from  the  topographic  surface  to  the  sphere  r  =  a.  This  can  be  done  by 
first  reducing  the  data  from  the  earth’s  surface  to  the  ellipsoid  or  the 
geoid  as  follows  (see  Rapp,  ibid.,  p.  15): 


Agj  =  Ag 


lig 

ar 


h  + 


-k 

2! 


-  ... 


(2.3) 


where  Ag  refers  to  the  earth’s  surface,  Ag^  refers  to  the  ellipsoid,  and 
h  can  be  considered  the  orthometric  height  of  the  surface  point.  As  a 
next  step  that  we  will  detail  later,  the  Ag^  will  be  continued  to  Aga 
referred  to  the  equatorial  sphere  (see  Figure  3). 

Operationally  the  gradients  in  (2.3)  can  be  obtained  from  an  existing 
high  degree  spherical  harmonic  expansion,  in  application  of  the  ideas  in 
Rapp  (ibid.,  p.  20).  The  correction  terms  in  (2.3)  are  then  evaluated 
numerically  and  applied  to  the  first  term  Ag  to  obtain  Ag^.  Since  in 
this  operational  solution  (2.3)  is  not  being  implemented  rigorously,  it  is 
appropriate  that  its  application  be  limited  to  the  continuation  of  Ag  to 
Agg  over  the  distance  h,  and  that  it  is  not  used  to  continue  Ag^  to  Ag, 
over  the  distance  (a  -  r^).  As  can  be  seen  in  Rapp  (ibid..  Table  1),  the 
magnitudes  of  the  correction  terms  to  transform  Ag  to  Ag^  are  on  the 
average  an  order  of  magnitude  smaller  than  the  terms  to  transform  Ag^ 
to  Aga. 


I 


For  the  continuation  of  Ag^  to  Aga  we  implicitly  use  also  a  Taylor 
series  expansion,  analogous  to  (2.3),  but  this  time  we  find  that  the 
Taylor  series  terras  can  be  evaluated  in  an  analytical  as  opposed  to  a 
numerical  fashion.  This  is  possible  because  the  deviation  between  the 
surfaces  of  the  ellipsoid  and  equatorial  sphere  is  such  a  simple  function 
of  latitude  that  the  Taylor  terms  themselves  can  be  expanded  in 
spherical  harmonics.  Therefore,  we  find  that  in  a  single  derivation  we 
can  express  the  unknown  potential  coefficients  in  terms  of  Ag^,  without 
any  need  to  continue  Ag^  to  Aga  is  a  separate  step.  The  starting  point 
of  the  derivations  is  (1.13)  applied  on  the  ellipsoid: 


e* 


-  *Tt 

sin^cos^  — ^  + 

rt** 


[eJa  Pa(sin^)  -  3  (1  -  sin**)]  Tj , 


(2.4! 


with  Te  obtained  by  applying  (2.1)  on  the  ellipsoid: 


T.  =  “2_  if-r  c«.wi.  X). 

n=2  in=-n  ' 


(2.5) 


Before  proceeding  to  the  detailed  derivation  of  C*,„  from  using  (2.4) 

and  (2.5),  we  should  note  that  (2.5)  implies  the  use  of  the  series  (2.1) 
below  its  sphere  of  guaranteed  convergence,  i.e.,  the  sphere  of  radius 
r  =  a.  However,  numerical  investigations  by  Jekeli  (1981)  using  a 

realistic  gravity  field  model  would  indicate  that  the  use  of  (2.5)  is 
justified  at  least  up  to  maximum  degrees  Nmn^  of  practical  interest,  e.g., 
around  300. 

Let  us  now  discuss  the  recovery  of  potential  coefficients 
occurring  in  (2.1),  from  an  inversion  of  (2.4)  with  (2.5)  given  Ag^.  In 
our  derivations  we  will  arrive  at  so-called  ellipsoidal  correction 
expressions.  A  check  of  our  derivations  is  provided  by  the  agreement 
of  certain  expressions  with  those  of  Pellinen  (1982)  who  also  deals  with 
ellipsoidal  corrections.  The  newness  of  our  results  over  those  of 
Pellinen  (ibid.)  consists  in  deriving  new  equations  (e.g.,  (2.33)-(2.36) 
giving  ellipsoidal  corrections  when  geocentric  latitudes,  not  geodetic 
latitudes,  are  used  as  coordinates)  and  in  using  a  different  derivational 
approach.  Our  derivations  use  elementary  Taylor  series  considerations 
and  naturally  lead  to  insightful  physical  interpretations  of  the  various 
formulas.  The  general  idea  of  the  solution  will  be  to  first  somehow 
transform  (2.4)  into  the  form: 


•  n 

‘"“EH”  nm?  ^nm  (?.  X) 
n=2  n=-n 


(2.6) 


where 

conventional  spherical  harmonic  coefficients;  these 
are  constants  independent  of  (?,  X);  the  subscript 
▼  is  used  to  denote  the  use  of  geocentric  latitude 
in  the  expansion. 

Ynm(^i  X)  conventional  spherical  harmonics;  the  use  of  con¬ 
ventional  as  opposed  to  fully-normalized  harmonics 
is  for  convenience  in  the  derivations. 


It  will  turn  out  after  the  transformation  of  (2.4)  to  the  form  (2.6)  that 
the  are  expressed  as  some  function  of  the  conventional  potential 

coefficients  Cn„  (unnormalized  counterpart  of  C*„),  symbolically; 


T 
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An  inversion  of  (2.7)  yields,  synbolically: 


Cnm  ~  ^*(DnmV) 


(2.8) 


Meanwhile,  the  D„aif  can  also  be  directly  obtained  by  an  inversion  of 
(2.6)  given  the  hgi  (see  Heiskanen  and  Moritz,  1967,  1-70): 


D„o*  =  ^  1  I  AgE(V.  »  Pn(sin  f)  Ar  ; 


^  I  J  Agf  (♦.  X)  ¥„,(♦.  X)  da  (2.9) 


(■  *  0;  k  =  l.l), 


Therefore,  the  combination  of  (2.8)  and  (2.9)  yields  the  desired  solution 
for  the  recovery  of  potential  coefficients  Cn,,  from  Ag^;  the 
full-normalization  of  yields  C*ai  (see  (2.19)  below)  which  is  the  final 
desired  quantity.  The  just  outlined  solution  involves  several  steps, 
which  we  will  now  describe. 

First,  let  us  collect  the  following  relations  where  terms  proportional 
to  e*  have  been  dropped  (e  is  the  earth’s  first  eccentricity,  e’«0.()067): 


rg  =  a  (1  - ^  sin  *  ♦)  +  0(e*) 


(f!l 


n+2  n+2  ,  . 

*  1  +  e*  sin^t 


(2.10) 


(2.11) 


(^1 


n+1 


»  1  +  e*  8in*4 


Pa (sin  ♦)  =  §  sin*#  -  j  (exact) 


«‘r 


kM  2a 


.  »  a  ' 

.v-‘N 

« , 
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777:^ 
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,s* 


(2.12) 

»  y*.' 

(2.13) 
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1967,  (2-118) 

(2.14) 

V " 

m 

(2-127)) 

(2.15) 

m 


Thirdt  let  us  re-write  the  various  terms  of  (2.4).  Usintf  (2.17)  and 
(2.11)  the  first  two  terms  of  (2.4)  become: 

'tIf  *  ^  ?  li  II  ^ 

n=2  ■=-n 

^  ^  IZ  IZ  Cn«  *in»?  Y„,(?.  X).  (2.20) 

n=2  ■=-n 

Note  that  for  the  purpose  of  (2.20)  we  used  ?  instead  of  4  in  the  (small) 
e^-term  of  (2.11).  To  re-write  the  third  term  of  (2.4)  we  substitute 
(2.17),  use  (2.10)  and  (2.12),  and  neglect  terms  involving  e*  to  get: 


.  .  _  _  *Tt  -e^kM  _  .  _  _  . 

-e*  sin  f  cos  ? - 1-  =  — j-  \  )  C„*  sin  f  cos  ?  ”  (▼,  X). 

n=2  ^=-0  (2.21) 


To  re-write  the  last  group  of  terms  in  (2.4),  we  first  use  (2.14),  (2.10), 
(2.13),  and  (2.15)  to  simplify  the  quantity  in  brackets  in  (2.4);  then  we 
use  (2.12)  and  (2.17);  in  all  manipulations  we  drop  terms  involving  e*; 

we  get: 


6J,  ^  Pa  (sin  ▼)  -  (1  -  sin»V)  Te  =  Si  (3  sin»?  -  2)Te 

(2.22) 


=  ^  ^  C„„  sin*?  Y„,.(?,  X)  - 

n=2  B=-n 


-  C„„  Y„„(?.  X).  (2.23) 

n=2  B=-n 

Fourth,  let  us  write  three  relations  that  have  always  been  found 
useful  in  ellipsoidal  correction  theory.  The  first  of  these  relations  is 
(see  Moritz,  1980,  39-46): 


(2.25) 


(2.26) 


(2.27) 


Note  that  the  left  hand  sides  of  (2.24)  and  (2.26)  occur  in  (2.20),  (2.21), 
and  (2.23),  and  this  is  why  (2.24)  and  (2.26)  are  useful  to  us.  Finally, 
the  third  useful  relation  is  a  simplifying  relation,  which  says  that  a 
summation  index  can  be  shifted;  for  example  (see  ibid.,  39-48): 


y  ~  ^nm  ^  ^n+a,ai  ^n+2,m  •  (2.28) 

n=2  n=0 


It  is  to  be  noted  that  the  use  of  relations  such  as  (2.28)  introduces  zero 
and  first  degree  (i.e.,  n=0  and  n=l)  harmonic  terms  in  the  formulas. 
However,  as  customary  in  gravimetric  analyses,  we  will  suppress  these 
terms  for  convenience  by  always  starting  all  summations  from  n=2.  The 
consideration  of  these  n=0,  1  terms  can  be  treated  as  a  separate 

problem,  such  as  what  was  done  in  Lelgemann  (1970),  with  numerical 
evaluations  in  Rapp  (1981a),  for  geoid  undulation  computations. 
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Fifth,  we  substitute  (2.20),  (2.21),  and  (2.23)  into  (2.4),  then  use  the 
relations  (2.24)  and  (2.26).  Summation  indices  are  then  shifted  as  in 
(2.28)  in  order  to  have  only  ¥„„(?,  X),  and  not  Y„+a,(iiffi  X)  or 
Y„-a  ^)i  occurring  in  the  expression.  All  coefficients  of  Ynn,(V,X), 
which  by  now  are  all  constants  independent  of  position  (?,  X)  are  then 
collected  to  yield  a  coefficient  which  must  then  be  exactly  equal  to  the 
D„a,g  of  (2.6).  Observing  further  that,  from  (2.25)  and  (2.27): 


1  ®  ^nm 

(2.29) 

1  *  “  1 

(2.30) 

=  (n+l)  , 

(2.31) 

We  finally  get  the  simplified  expression: 


Dn-f  =  “  (n-l)C„*  +  ^  {(n^-n-^2)  C„- 


+  [(n*+n+l)^„a,-3]  +  (n*+3n+4)  rn+»,<ii  C„+a,i»} 


(2.32) 


which  is  in  the  form  (2.7). 

The  last  equation  agrees  with  Pellinen  (1982,  eq.  14),  after  putting 
m=e*/2  (cf.  Moritz,  1980,  39-12).  As  in  Pellinen  (1982),  underlining  a 
term  will  meeui  that  the  term  will  be  omitted  for  n<4  (note  in  (2.32)  that 
if  n=3  for  example,  the  underlined  term  will  involve  which  is  to  be 

suppressed  since  we  are  suppressing  all  harmonics  with  n<2). 

Given  the  values  Ag^  on  the  ellipsoid,  the  form  (2.6)  can  actually  be 
inverted  to  yield  the  constants  D„„,g  directly.  In  fact,  the  inverted 
expression  is  (2.9).  Therefore,  we  can  assume  that  the  on  the 

left-side  of  (2.32)  is  known,  and  hence  we  are  now  in  a  position  of 
being  able  to  recover  the  potential  coefficients  occurring  in  (2.17), 
simply  by  inverting  (2.32).  To  do  the  inversion,  in  the  first  place  we 
set  e^=0  in  (2.32)  to  get  the  sphericed  approximation  solution: 


^nmW  ~ 


kM(n-l) 


(2.33) 


where  the  superscript  denotes  spherical  approximation  and  the 

subscript  denotes  the  use  of  geocentric  latitude.  Knowing  the 
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these  values  are  used  in  the  e^  terms  of  (2.32),  and  the  equation  then 
inverted  to  yield  the  Cn„.  In  addition  uaing  (2.27),  the  solution  can  be 
finally  written  in  the  form: 


“  ^nni¥ 


where 


^^nmV  ”  ®*(PnmV  ^n— "*■  ^nniV  ^n+2,BiV)t 


(2.34) 


(2.35) 


and  with  k  =  Iml 


(n-k-1) (n- 


""  2(n-l)(2n-3)(2n-l) 


-2n*>2n-«k^-4n»-(-2nk^->-9n»+2k»-»-lln-8 
qnm?  -  2(n-l)(2n+3)(2n-l) 


-(n*+3n+4)  (n->-k-*-2)  (n-t-k-*-] 
2(n-l)(2n+5)(2n+3) 


(2.36) 


The  above  equations  agree  with  the  ones  recently  given  by  Rapp  (1984). 

In  some  cases  we  may  be  handling  the  Ag^-data  using  geodetic 
latitudes  ^  instead  of  geocentric  latitudes  7.  In  such  cases  we  have, 
instead  of  (2.6),  the  expansion: 


4gE  =  Dnm^  ^nm  (♦.  X) 

n=2  m=-n 


(2.37) 


We  now  want  to  know  how  the  solution  consisting  of  (2.33)-(2.36)  will 
change.  First,  let  us  relate  and  Dn„f  by  transforming  (2.6)  to  the 

form  (2.37).  Omitting  terms  involving  e*,  we  have  from  Rapp  (1981a, 
p.  2)  that: 


Y„m(^,  X)  =  ¥„„(♦,  X)  + 


(▼  -  ♦)  ®  -e*  sin  W  cos  ?  . 


(▼  -  ♦)  : 


(2.38) 


(2.39) 


We  substitute  (2.38)  and  (2.39)  into  (2.6),  then  use  the  relation  (2.24). 
Analogous  to  the  derivation  of  (2.32),  summation  indices  are  then  shifted 
in  order  to  have  only  Y„„,  and  not  Y„-a,m  or  Y„_2,m,  occurring  in  the 
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expression.  All  coefficients  of  X),  which  by  now  are  all  constants 

independent  of  (^i  X)i  are  then  collected  to  yield  a  coefficient  which 
must  be  exactly  equal  to  the  coefficient  of  (2.37).  We  get: 

^nin4  "  ®  (3n— 2,<n  ^n— 2,ni  l>nm  lln  +  2,m)  (2.40) 

which  is  the  desired  relation  between  and  D„„f.  Substituting 

(2.32)  into  (2.40),  using  the  approximation  (see  (2.33)) 


Dn«  *  ^  (n-1)  C„„ 


(2.41) 


in  the  (small)  terms  involving  e*,  and  simplifying  we  get: 
Dn«.4  =  ^  (n-1)  C„„  +  ^  ((3n»-lln-H4)  c.„-,  ,  C„-,  - 
+  [(n^+n+l)  -  2(n-l)  b„„  -  s]  C„^ 

-  (n*+5n+2)  rn+2,m  C„+,,4  . 


(2.42) 


Equation  (2.42)  can  now  be  inverted,  in  exactly  the  same  way  as  (2.32) 
was  inverted,  to  give  the  following  solution  based  on  geodetic  latitudes: 


C„m4  = 


kM(n-l) 


^nm  -  ^ntn4 

where 

“  ®*(Pn«n4  ('n-2,n<4  *1nm4  ^nint  >^nm4  ^n+2,in4), 

and  with  k  =  Iml : 


(2.43) 


(2.44; 


Pniii4 


-  -(3n^-lln-H4)(n-k-l)(n-k) 
2(n-l)(2n-3)(2n-l) 


*lnm4  " 


-2n*+2n'k*-2n*-4nk*+9n*+8k*+9n-8 

2(n-l)(2n+3)(2n-l) 


(2.45) 


_  (n^-*-5n-»-2)  (n-»-k>2)  (n-*-k-*-I ) 
■  2(n-l)(2n-^5)(2n+3) 


‘.Vi 

>>> 


;vS 


S3 


O  •s.* 


*  _• 


‘v 

>?/ 

8::^ 


IT: 
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The  above  equations  agree  with  Pellinen  (1982,  (16)-(18)). 

In  summary,  we  are  given  the  values  of  gravity  anomalies  Ag^  on 

the  surface  of  the  ellipsoid.  The  Ag^  are  the  analytically  continued 
version  of  the  actually  available  Ag  on  the  earth’s  surface.  The 
anomalies  on  the  ellipsoid  can  be  formally  expanded  into  a  series  (2.6)  or 
(2.37)  of  spherical  harmonics,  with  coefficients  or  obtained 

through  (2.9).  The  coefficients  Dnm^  or  can  then  be  used  in 

(2.43)-(2.45)  or  (2.33)-(2.36),  respectively,  to  generate  the  conventional 
potential  coefficients  (aee  (2.17)).  Finally,  the  C„„  can  be  fully 

normalized  by  (2.19)  to  yield  the  C*„  occurring  in  (2.1),  completing  the 
recovery  of  potential  coefficients. 

In  our  derivations  we  have  always  omitted  terms  involving  the 
fourth  and  higher  powers  of  the  earth’s  first  eccentricity  e.  This  does 
not  mean,  however,  that  the  relative  error  introduced  by  omission  of 
terms  is  only  0(e*)  every  time.  Specifically,  in  (2.11)  and  (2.12)  (see 
comments  under  (2.15))  the  relative  error  depends  also  on  the  harmonic 
degree  n  and  on  the  latitude  4,  in  addition  to  depending  on  the  omitted 
powers  of  e.  In  the  case  of  (2.12),  the  omission  of  terms  was  not 
critical,  and  in  fact,  only  the  first  term  of  (2.12)  was  eventually  used; 
the  reason  is  that  the  use  of  (2.12)  was  to  write  (2.21)  and  (2.23),  which 
were  equations  that  already  involved  only  small  terms  of  0(e*).  In  the 
case  of  (2.11),  however,  the  omission  of  terms  was  critical;  the  reason  is 
that  (2.11)  was  used  to  write  (2.20),  which  was  an  equation  for  the 
(principal)  first  two  terms  of  (2.4).  In  fact,  the  dominating  source  of 
inaccuracy  in  our  whole  derivations  was  the  omission  of  terms  in  (2.11), 
directly  causing  significant  relative  error  to  (2.20)  (see  error  estimates 
later  this  section). 

For  higher  accuracies  we  need  to  consider  the  next  omitted  term  in 
(2.11),  namely,  the  term  (2.16).  This  can  be  done  as  follows.  First, 
there  will  be  the  following  additional  term  to  the  right-side  of  (2.20): 


^  Y2  IZ  ("■^2)(n+l)(n-l)  sin*V  (▼.  \)  . 
n=2  in=-n 


By  recursive  application  of  (2.26)  we  have: 


(2.46) 


sin*?  Y„„(f,  X)  =  + 


(2.47) 


where 


~  ®niii  *n+a,m 


K„. 

Lnoi  “  *n(n(^n«  Pn+t,m) 

Mnm  =  “nm  rn+a.m  ■*■  ■*■  *n-ai<"  '>'"1" 

Num  “  TffifliC^n^3|M 
Pnm  =  7n-2,m  7nm  • 


(2.48) 


Equation  (2.47)  with  (2.48)  ia  also  given  in  Pellinen  (1982,  (9)). 

Analogous  to  the  construction  of  D„„g  in  (2.32),  we  substitute  (2.47)  into 
(2.46).  Summation  indices  are  then  shifted  as  in  (2.28)  in  order  to  have 
only  Y„„,  and  not  (Y„+4,„,  Y„+,,«,  Y„.,,»,  Y„,*„)  occurring  in  the 


expression, 
coefficient  D„a,9 


All  coefficients  of 

,1  • 


Y„«(?,  X)  are 


then  collected  to  yield  a 


e^kM 


"nnif  -  8a* 


(2.49) 


where 

SnmV  ”  (n~2)  (n~3)  (n~5)  K„— 4^,,  C„— 4^,1,  + 
+  n(n-l)(n-3)  L„-3,ni  Cn-a.m 
+  (n+2)(n+l)(n-l)  Mnm  ^nm 
+  (n+4)  (n+3)  (n+1)  N„-f2,ni  Cn+a,m 
+  (n+6) (n+5) (n+3)  Pn+4,m  ^n+4,w  • 


(2.50) 


(2.51) 


A  more  complete  version  of  (2.32)  can  now  be  written  as  follows: 


„  kM  ,  .  e^kM  „ 

Dnm?  =  P"  +  gg* 


e*kM 


(2.52) 


where  R„«s  denotes  the  quantity  in  braces  in  (2.32).  Equation  (2.52) 
can  then  be  inverted  to  yield  the  unknown  C„„.  The  inversion  can  be 
done  by  first  writing  (2.52)  in  the  form: 


kM(n-l) 


OnmV 


2(n-l) 


_  _ e: 


8(n-l) 


5nmV 


(2.53) 
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From  (2.53)  we  have  the  "zero  order"  solution: 
_ a* 


kM(n-l) 


(2.54) 


Equation  (2.54)  is  the  same  as  the  spherical  approximation  solution 
(2.33),  and  corresponds  to  the  case  when  only  the  first  term  in  (2.11) 
(therefore,  only  the  first  term  in  (2.20))  is  used.  Usin^  (2.54)  in  (2.53) 
we  obtain  the  "first-order"  solution: 


'nm  =  _  2(n-l) 


,_i  \  ®n<iiV 


(2.55) 


Equation  (2.55)  is  the  same  as  the  solution  (2.34),  and  corresponds  to 
the  case  when  the  first  and  second  terms  of  (2.11)  (therefore,  also  the 
first  and  second  terms  in  (2.20)),  plus  the  correction  terms  of  O(e^) 
tfiven  by  (2.21)  and  (2.23),  are  accounted  for  in  the  solution.  Finally, 
using  (2.55)  into  (2.53)  yields  the  "second-order"  solution: 


.(»)  _ 


=  C 


2(n-l)  8(n-l) 


sLv 


(2.56) 


where  the  superscript  "1"  means  that  C‘_4,„  C*-*,,,,  C‘„,  C‘+*  _ 

from  (2.55)  are  used  to  obtain  and  The  solution  (2.56) 

corresponds  to  the  case  when  the  terms  of  (2.11)  and  (2.16)  (therefore, 
also  the  terms  of  (2.20)  and  (2.46))  are  accounted  for  in  the  solution 
along  with  the  terms  of  (2.21)  and  (2.23). 


A  physical  interpretation  can  be  given  to  the  zero-,first-,  and 
second-order  solution  for  the  recovery  of  potential  coefficients  from 
Age-data.  This  is  seen  from  first  noting  that  we  can  write  the  Taylor 
expansion: 


d«E  =  dg. 


I  (rf-a)*  + 
'  r=a 


_1  «Mg| 


3!  «r= 


'  r=a 


(rg-a)* 


(2.57) 


where  we  define: 


dg  ■  -  -  T 

®  *r  r 


^  E  E  (?) 

n=2  m=-n 


n+2 


’(n-1)  Cn„  Y„„(V,  X) 


(2.58) 


J*"'  w"*  •*“  ■**>  '  »*•  .  ”  . 


% 


12 


V.'.- 


ii 


•,*  s' 


rf 


30 


(rg  -  a)  =  "  -j-  sin*4  +  0(e*)  ; 

dgE  =  dg(rE,  ?,  X)  ;  (dg  on  the  ellipsoid) 

dga  =  dg(a,  ?,  X)  .  (dg  on  the  equatorial  sphere) 


Equation  (2.20)  with  the  additional  term  (2.46)i  which  is  a  starting  point 
of  our  derivation,  is  in  the  general  form  of  the  Taylor  expansion  (2.57) 
truncated  after  the  third  term.  The  zero-order  solution  in  effect  uses 
only  the  first  term  in  (2.57),  thereby  setting  dga^dg^;  this  means  that 
in  this  solution  a  given  boundary  value  Ag^  at  position  (re,  ▼>  X)  is 
formally  considered  to  be  on  the  equatorial  sphere  at  position  (a,  ?,  X). 
The  first-order  solution  on  the  other  hand,  uses  the  first  two  terms  of 
(2.57),  thereby  using  the  first  radial  derivative  of  the  dg-field  to 
distinguish  between  dga  and  dg^.  Finally,  the  second-order  solution 
uses  the  first  three  terms  of  (2.57),  thereby  using  the  first  and  second 
radial  derivatives  of  the  dg-field  to  distinguish  between  dga  and  dg^. 

To  examine  numerically  the  effect  of  omission  of  terms,  we  write 
(2.11)  more  completely  as: 


^ 


"2  sin*4  + 


j+1)  e' 


—T  sin*4  + 
4 


-g  sin*4  + 


(2.62) 


With  e‘/2  0.003  -  f  (the  earth’s  flattening),  the  maximal  orders  of 

magnitude  of  the  terms  in  (2.62)  are  obtained  at  the  poles  with  sin4  =  1: 


(2.63) 


0(e*-tenn)  =  nf  =  3  x  10  ’  n 


0(e*-tenn)  =  =  5  x  10  ‘  n* 


0(e*-tenD)  =  “o7~  =  5  x  10“* 


The  omission  of  terms  estimated  in  (2.63)  directly  causes  a  relative  error 
in  our  derivations.  Since  this  source  of  error  is  dominating,  we  can 


(2.59) 

P 

(2.60) 

A"!** 

(2.61) 

►  “  •  ’ 
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alao  interpret  (2.63)  as  giving  maximal  estimates  of  the  total  relative 
error  itself,  caused  by  the  omission  of  certain  powers  of  e  in  the 
recovery  of  potential  coefficients  €*„  from  Age-data.  Such  estimates, 
which  agree  with  those  of  Pellinen  (1982),  are  eveduated  in  Table  1  for 
various  values  of  n  and  various  powers  of  omitted  e-terms. 

Table  1  says  that  for  n=200,  for  example,  the  omission  of  e’  and 
higher  powers  (i.e.,  the  use  of  only  the  first  term  in  (2.34)j  (2.44),  or 
(2.53))  causes  a  relative  computational  error  of  60%  in  Cfoo^mi 
percentage  error  decreases  to  20X  and  to  4X  when  terms  proportional  to 
e*  and  e*,  respectively,  are  considered  in  the  formulas. 


Table  1.  Relative  Error  Caused  by  the  Omission  of  Terms  Proportional 
to  Certain  Powers  of  e,  in  the  Recovery  of  Potential 
Coefficients  cHai  from  Anomaly  Data  Age  on  the  Ellipsoid. 


MiniauB 
Power  of  e 
in  Oaitted 
Terms 


Relative 
Error  cr 

Value  of  Relative  Error 
for 

f^OiACS^l 

L  0(C*,„)  J 

n=30 

n=100 

n=200 

n=300 

3  X  10"’  n 

0.09 

0.03 

0.60 

0.90 

5  X  10-‘  n» 

0.0045 

0.05 

0.20 

0.45 

5  X  10-»  n=» 

0.0001 

0.005 

0.04 

0.13 

The  relative  errors  cr  given  in  Table  1  affect  the  gravimetric 
quantities  computed  from  the  C*,,.  We  can  look  at  short  and  long 
wavelength  effects  separately  by  looking  at  effects  on  computed  gravity 
anomedies  Ag  and  height  anomalies.  The  global  RMS  effects  when  using 
an  expansion  to  degree  and  order  N^ax  can  be  expressed  as  the  square 
root  of  (in  a  sphere  of  radius  r=a): 


N,., 

=  (^1*  [0(«C*„ 


r»(f)  =  a*  (2n+l)  [0(AC*  )]»  ; 


(2.64) 


(2.65) 


ck  is  the  relative  error  from  Table  1  or  (2.63)  . 


(2.67) 


From  Kaula’s  rule  of  thumb,  we  have  the  following  order  of  magnitude  of 
potential  coefficients: 


0(cL)  =  ^ 


(2.68) 


Let  us  further  use  values  from  the  Geodetic  Reference  System  1980: 


kM  =  3986005  x  10*  m^  s"* 
a  =  6378137  m  . 


(2.69) 

(2.70) 


Using  (2.64)-(2.70),  values  of  <r(4g)  and  <r(()  were  computed  and  shown 
in  Table  2,  for  various  maximum  degrees  of  expansion  N*,,.  For 
Nm«x=200,  for  example,  the  effect  of  omitting  the  (e*,  e*,  e*)  and  higher 
powers  would  be  to  cause  an  error  of  about  (5.9,  1.4,  0.2)  mgals  in 
grt  dty  anomalies  and  an  error  of  about  (0.62,  0.06,  0.01)  meter  in  height 
anomalies,  on  a  sphere  of  radius  r=a. 

Table  2.  Effect  of  Potential  Coefficient  Errors  on  Computed 

Gravity  Anomalies  4g  and  Height  Anomalies  (. 

Units:  <r(4g)  in  milligals;  <t(0  in  meters. 


Minimum 
Power  of 
e  in 
Omitted 
Terms 


Relative 
Error  Cfi 
of  Potential 
Coefficients 
(from  Table  1) 


Value  of 

[<T(Ag)l 

^e(f) 

from 

[(2.64)1 

U2.65)J 

5  X  10"*  n 


N„,x=30  N„„=100  N„„=200  N„„.=300 


0.85  mgals 

2.91 

5.85 

8.79 

0.49  m 

0.57 

0.62 

0.64 

0.03 

0.35 

1.39 

3.12 

0.01 

0.03 

0.06 

0.10 

5  X  10-*  n» 


0.00 

0.00 


0.03 

0.00 


0.77 

0.02 


It  is  not  surprising  that  relatively  large  errors  are  associated  with 
the  zero-order  solution  (see  the  first  row  of  values  in  Tables  1  and  2): 
the  reason  is  that  the  true  gravity  anomaly  4g,(a,  J,  X)  on  the 
equatorial  sphere  and  Ag^lrE,  V,  X)  on  the  ellipsoid  can  easily  differ 
significantly,  so  that  equating  ig,=AgE  (strictly,  equating  dgaCdg^),  as 
the  zero-order  solution  does,  incurs  large  errors.  The  consideration  of 
the  e’-terms,  as  in  the  first-order  solution,  means  the  consideration  of 
the  first  radial  derivative  of  the  gravity  anomaly  field  to  relate  4g,  and 
resulting  in  much  smaller  errors  (see  second  row  of  values  in 
Tables  1  and  2)  than  those  of  the  zero-order  solution.  Finally,  a 
consideration  of  e*  terms,  as  in  the  second-order  solution,  means  the 
additional  consideration  of  the  second  radial  derivative  of  the  field  to 
relate  Ag,  and  Agg,  resulting  in  even  smedler  computational  errors  (third 
row  of  Tables  1  and  2).  As  seen  from  Tables  1  and  2,  the  higher  the 
harmonic  degrees  (n)  of  interest,  the  higher  the  order  of  radial 
derivatives  of  the  gravitational  field  that  need  to  be  considered  in  the 
ellipsoidal  correction  theory. 

Summarizing,  the  determination  of  potential  coefficients  requires  the 
surface  gravity  anomalies  to  be  analytically  continued  to  some  sphere. 
Rather  than  continuing  to  a  bounding  sphere  as  in  Rapp  (1984),  we 
recommend  continuing  to  the  equatorial  sphere  where  the  standard  form 
of  spherical  harmonic  representation  (2.1)  refers.  For  operational 
implementation  we  recommend  carrying  out  the  analytical  continuation  in 
two  steps.  The  first  step  would  be  to  continue  the  surface  data  to  the 
ellipsoid.  This  can  be  done  (ibid.)  by  numerically  computing  the 
correction  terms  of  the  Taylor  series  using  surface  heights  and  vertical 
gradients  from  an  existing  spherical  harmonic  expansion.  The  second 
step  would  be  to  analytically  continue  from  the  ellipsoid  to  the 
equatorial  sphere.  This  step  can  be  carried  out  also  by  Taylor  series, 
but  with  the  Taylor  terms  applied  analytically  through  the  ellipsoidal 
corrections  of  Pellinen  (1982).  At  maximum  degrees  of  expansions  around 
n  =  300,  up  to  second  order  Taylor  terms  need  to  be  applied,  and  this 
can  be  done  through  (2.56).  Tables  1  and  2  indicate  the  number  of 
Taylor  terms  required  for  given  maximum  degrees  of  expansion  and 
desired  accuracy  of  continuation  from  ellipsoid  to  equatorial  sphere. 
The  effects  of  the  Taylor  terms  for  the  continuation  from  the 
topography  to  the  ellipsoid  are  on  the  average  an  order  of  magnitude 
less  than  the  values  in  Tables  1  and  2  (see  Rapp,  ibid..  Table  1). 


3.  FAIIILIASIZATION  WITH  THE  SPATIAL  DISTUHBANCB  VECTOR  SIGNAL 


Given  potential  coefficients  the  anomalous  potential  T(r,  X)  is 

expressed  in  spherical  harmonics  by  (2.1).  Then  with  (2.1)  as 
fundamental  equation,  other  gravimetric  quantities  of  interest  can  be 
derived  in  spherical  harmonics  using  (2.1)  and  its  first  and  second 
order  derivatives.  As  an  example,  Rapp  (1982)  gives  a  collection  of  such 
derivations  for  common  gravimetric  quantities,  specifically  the  height 
anomaly,  gravity  anomaly,  three  components  of  the  gravity  disturbance 
vector,  and  two  components  of  the  deflection  of  the  vertical.  On  the 
other  hand,  given  the  gravity  anomalies  Ag^  on  the  ellipsoid,  the 
anomalous  potential  T(r,  V,  X)  is  expressed  in  an  integral  formula  by 
(3.14)  below.  Now  using  (3.14)  and  its  first  and  second  order 
derivatives,  integral  formulas  for  other  gravimetric  quantities  can  also 
be  derived.  Examples  of  such  derivations  are  found  in  Heiskanen  and 
Moritz  (1967,  Sec.  6>4)  for  the  components  of  both  the  gravity 
disturbance  vector  and  the  deflection  of  the  vertical. 

Both  parameterizations  from  0*^  and  Ag^  have  their  respective 
theoretical  and  practical  problems.  The  theoretical  problem  with  the 
parameterization  is  the  divergence  of  downward  continued  spherical 
harmonic  series.  The  total  (i.e.,  infinite)  spherical  harmonic  series  (2.1) 
representing  the  anomalous  potential  is  not  guaranteed  to  converge,  and 
in  fact  will  probably  diverge,  when  evaluated  below  the  smallest  sphere, 
called  the  minimum  sphere,  bounding  the  earth.  The  errors  arising  from 
this  total  series  divergence  generally  grow  with  the  depth  of 
penetration  into  the  minimum  sphere  and  the  power  (expressed  by 
degree  variances)  of  the  high  frequency  variations  in  the  anomalous 
field.  For  depths  of  penetration  down  to  the  earth's  surface  and  for 
the  degree  variance  decay  of  the  earth's  anomalous  potential,  Jekeli 
(1981)  shows  completely  negligible  errors  coming  from  total-series 
divergence  when  using  a  spherical  harmonic  representation  of  the 
earth’s  anomalous  field  to  maximum  degree  N„„=300. 

Of  great  concern  also  would  be  the  practical  problems,  associated 
with  the  parameterization.  The  slow  decay  of  the  cj,,,  with  n 

necessitates  the  handling  of  a  very  large  number  of  coef^cients  in 
order  to  represent  all  significant  variations  in  the  anomalous  field.  For 
example,  a  S'xS'  resolution  which  may  be  needed  in  certain  computations 
necu*  the  earth’s  surfeu:e  would  require  approximately  Nn,a„=2000  or  a 
total  of  (2001)^^4  million  coefficients.  The  use  of  such  high  degree 
spherical  harmonic  representation  will  require  a  prohibitive  amount  of 
computer  storage  and  time.  Also,  such  high  degree  representation  will 
require  some  sort  of  re- parameterization  from  the  available  Ag  to  the 
and  such  re-parameterization  will  involve  difficult  practical 


problems  related  to  computer  requirements,  quadratures,  and 
implementation  of  a  relation  between  C*hi  and  that  accounts  for  the 
elUpticity  and  topography  of  the  earth  (Rapp,  1984).  Finally,  there  is 
the  problem  of  updating  a  spherical  harmonic  representation  to 
incorporate  updated  gravity  anomaly  data.  Such  updates  are 
inconvenient  as  they  everytime  require  a  new  solution  for  the  entire 
set  of  parameters.  Today,  taking  altogether  the  theoretical  and  practical 
problems  of  the  C*,„-parameterization  it  would  appear  that  a  feasible 
upper  limit  to  use  would  be  around  N„ax=360. 

At  higher  frequencies  (harmonic  degrees  greater  than  n=360  or 
spatial  resolution  greater  than  30'x30')  the  direct  parameterization  from 
dgE  via  an  integral  formula  based  on  (3.14)  below  would  be  appropriate. 
Since  the  high  frequency  variations  are  attributable  to  Ag^-data  close  to 
the  computation  point,  the  integral  formula  can  be  truncated  to  within  a 
data  cap  of  limited  radius  around  the  computation  point.  This  minimizes 
the  practical  problem  of  integral  formula  representations,  that  of  having 
to  integrate  data  over  a  large  cap. 

Since  the  spherical  harmonic  and  integral  representations 
complement  each  other's  deficiencies,  the  ideal  would  be  to  combine  the 
two  representations.  Such  combination  is  performed  under  the  heading 
of  "truncation  theory"  as  originated  by  Molodensky  (see  Molodensky,  et 
al.,  1962),  or  more  recently,  under  the  heading  of  "least  squares 
spectral  combination  theory"  (Moritz,  1975;  Sjoberg,  1981;  Wenzel,  1982). 
We  will  discuss  the  combination  by  truncation  theory  in  this  chapter, 
with  special  emphasis  on  the  computation  of  the  gravity  disturbance 
vector  components  at  altitude. 


3.1  Formulas  for  Spherical  Harmonic  and  Integral  RepreseniationB 
From  (2.1)  we  can  express  the  anomalous  potential  as: 

^ “  r  E  IT"  ‘)  ■ 

n=2  B=-n 


Using  (3.1)  in  the  usual  definitions  of  the  radial  and  horizontal  gravity 
disturbance  vector  components  and  the  definition  of  the  gravity  anomaly 
in  spherical  approximation,  we  can  collect  the  following  equations: 

Y~  y~  (n+1)  C„*  ¥„„(▼,  \) 

n=2  B=-n 


-«T  kM 
r  ~  «r  ~ 


(3.2) 
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. . ,  surface  hamonic  of  the  gravity  anomalies  on  the  geo¬ 
centric  sphere  of  radius  r^,  under  a  (▼,  X)-mapping  of 
these  anomalies  onto  a  unit  sphere. 


A  second  set  of  spectral  forms  can  be  written  analogous  to  (3.6)-(3.9)  in 
anticipation  of  the  use,  under  the  integral  representation,  of  the 
mapping  of  points  and  functions  discussed  in  Section  1.2  (see  Figure  2); 
we  have: 


n=2 

n=2 


=  E  1 

It?!”*' 

1  1  r 

•Sx 

n=2 

COS?  »X 

where  now: 


(3.10) 


(3.11) 


(3.12) 


j  ri=R+h  (3.13) 

R  ...  mean  earth  radiiis  (6371  km) 

'  h  . . .  height  of  the  computation  point  above  the  ellipsoid  or  the 

geoid 

I  •  •  jjth  surface  harmonic  of  the  gravity  anomalies  on  the 

ellipsoid  or  geoid,  under  a  (4,  X)-mapping  of  these  anomalies 
onto  a  unit  sphere.  I.e.,  geodetic  latitudes  ^  are  to  be  used 
in  the  defining  equation  (3.42)  below. 

f  .  Our  practical  use  for  the  type  of  equations  (3.10)-(3.12),  as  well  as  for 

*  their  versions  where  truncation  coefficients  have  been  introduced 

''  (details  next  sections),  will  only  be  for  degree  variance  propagations. 

This  is  in  view  of  the  fact  that  degree  variances  are  usually  quoted  for 
the  AgEn  of  (3.10)-(3.12)  and  not  for  the  &g„  of  (3.6)-(3.8).  For  actual 
\  computation  of  contributions  from  potential  coefficients,  however,  we  will 

I  use  the  type  of  equations  (3.6)-(3.8)  along  with  their  versions  where 

^  truncation  coefficients  have  been  introduced. 


i 

•• 


For  use  under  the  integral  formula  representation,  we  need  the 
space  domain  equivalents  of  (3.10)-(3.12).  These  space  domain 

equivalents  will  express  the  gravimetric  quantities  in  terms  of  The 

space  domain  expressions  for  (3.6)-(3.8)  follow  immedutely  from  those  of 
13*10)*{3*12)  b y  just  c hanging  to  r ^ ,  r ^  r ,  ^  to  ^ ,  and  A g ^  to  Ag  in 
the  formulas:  however,  such  expressions  will  be  essentially  of  no 
practical  value  since  we  do  not  have  the  required  detailed  Ag  on  a 
geocentric  sphere  of  radius  r^.  We  have-  the  following  space  domain 
equivalents  of  (3.10)-(3.12)  from  Heiskanen  and  Moritz  (1967,  Section 
6-4): 


T(r,  =  4“  J  J  S(ri,  V)  Agf  d<r 

<J 


(3.14) 


fir(r.  f,  X) 


4gE 


da 


(3.15) 


^X 


(r. 


cos  a 
sin  a 


da  , 


(3.16) 


where,  assuming: 

(4,  X):  geodetic  corrdinates  of  the  computation  point  P 

(4',  X'):  geodetic  coordinates  of  the  moving  integration  point  P', 

we  have: 

da  =  cos4'd4'dX' 

(integration  element) 

cow^  -  sin4sin4'  *  cos4co84'cos(X'-X) 

(angular  distance  between  P  and  P')  (3.17) 


tan  a 


cos  ♦'  3in(\'-X) 

cos  4  sin  4'  -  sin  4  cos  ♦'  cos(X'-X) 


(azimuth  from  P  to  P'). 


The  integral  kernels  in  (3.14)-(3.16)  can  be  obtained  in  the  space  domain 
using  (ibid.,  p.  235);  see  also  Paul,  1983,  (1);  Shepperd,  1982,  p.  102): 
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A  simple  check  is  provided  by  the  fact  that  if  we  take  the  negative  of 
the  derivative  of  (3.25)  with  respect  to  ri  we  arrive  at  (3.26)  as  it 
should  be.  The  spectral  form  of  the  function  (3.20)  which  forma  part  of 
the  kernel  for  the  horizontal  disturbance  components  (3.16)  can  be 
obtained  by  differentiating  (3.25)  directly,  and  verifying  from  the 
definitions  in  Heiskanen  and  Moritz  (1967,  p.  22)  that: 


7^  P„(y)  =  -Pni(y)  . 


(3.27) 


where  the  Pni(y)  are  the  associated  Legendre  functions  of  the  first 
order  (m=l).  We  arrive  at: 


E  (t?) 

n=2 


(3.28) 


Equation  (3.28)  has  a  similar  structure  to  (3.25)  and  (3.26),  except  that 
in  (3.28)  the  P„i(y)  and  not  the  Pn(y)  is  used  -  this  is  a  consequence 
of  the  non-isotropy  of  the  integral  kernel  in  (3.16). 


3.2  Spectral  Information  Content 

This  section  is  a  study  of  the  spectral  information  content  of  the 
radial  and  horizontal  gravity  disturbance  vector  components  at  various 
altitudes.  The  global  RMS  (root  mean  square)  value  of  the  radial 
disturbance  in  (3.11)  can  be  expressed  analogously  to  Heiskanen  and 
Moritz  (1967,  pp.  260-261)  as  the  square  root  of: 


=E  ■ 


(3.29) 


where  the  are  the  degree  variances  of  the  gravity  anomalies  Agg.  Of 
interest  also  is  the  RMS  information  beyond  a  given  harmonic  degree  N, 
expressed  from  (3.29)  as  the  square  root  of: 


(3.30) 


n=N+l 


parry  y>'*.v  \u  w,wj8  w.i'  ■  ^  ^  Jt  *3*  V*y*i^Jl  * J  *.f J  M  ^  M  M  PJ  nF  rj  rj  rip  v  v  vl""^  n;^vi'  ^  wyjr y ;  tv  ^  w  a?*  v 
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The  relative  omission  error  incurred  by  truncating  the  signal 
representation  to  a  maximum  degree  N  can  be  expressed  in  percentage 
as: 


esCfir)  =  X  100\  , 


(3.31) 


where  (TN(^r)  ^>^<1  <r(£r)  are  the  square  roots  of  the  quantities  <rf|^(6r) 
and  <i^{6r)t  respectively.  Quantities  analogous  to  (3.29)-(3.31)  can  be 
written  for  the  total  horizontal  disturbance  defined  as: 


6„»  =  V  +  6x*  . 


(3.32) 


where  6-^  and  6\  are  given  by  (3.12).  Applying  ibid.,  p.  262,  the  global 
RMS  value  of  the  horizontal  disturbance  in  (3.32)  can  be  expressed  as 
the  square  root  of: 


cfcSJ 


S'/- 


. .  ■ 


(3.33) 


n=2 


then  we  have  analogous  to  (3.30)  and  (3.31): 


"Z _  (n-l)»  IrJ  ’ 

n=N+l 


(3.34) 


(3.35) 


We  numerically  evaluated  the  above  expressions  for  various  altitudes 
h  above  the  ellipsoid,  using  R=6371  lim,  r,=R+h,  and  the  so-called 
two-component  degree  variance  model  (Rapp,  1979;  unit:  mgal’): 


C„  =  3.405  (0.998006)0+2  +  (0.914232)0+2  , 


n  k  3 


(3.36) 


Cj  =  7.5  ngal*  . 


••lyi 


>r 

E. 


I 
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The  summations  were  carried  out  to  n=5400,  corresponding  to  a  space 
domain  solution  resolution  of  2'  half-wavelength  (180*/5400=2' ). 

A  presentation  of  the  results  is  shown  in  Figures  4  and  5.  Since 
practically  identical  graphs  were  obtained  for  the  radial  as  for  the 
horizontal  component,  only  graphs  for  the  radial  component  are  shown. 
The  only  notable  difference  between  radial  and  horizontal  results  was 
that  ff(6r)  always  exceeded  <r(6n)  by  about  1.5  mgals  for  all  altitudes 
tested.  Figure  4  actually  shows  the  quantity  7N(5r)  from  (3.30),  or 
basically  also  the  quantity  from  (3.34).  The  parameters  used  were 

altitudes  h=S,  10,  20,  30,  50,  100,  200,  500  km,  and  various  harmonic 

degrees  defining  the  lower  limit  of  the  summation  in  (3.30).  For 

example,  Figure  4  says  that  the  use  of  a  radial  or  horizontal  disturbance 
representation  to  harmonic  degree  n=360  leaves  about  (16.2,  10.2,  4.5) 
mgals  to  be  resolved  at  (5,  10,  20)  km  altitude.  For  altitudes  of  100  km 

and  higher,  however,  such  representation  to  n=360  is  seen  to  be  essen¬ 

tially  complete  in  representing  all  the  information  content  of  the  field. 
A  representation  to  n=360  may  be  interpreted  as  .being  a  spherical  har¬ 
monic  representation  with  potential  coefficients  C*,,  complete  to  degree 

and  order  360,  or  it  may  be  interpreted  as  an  integral  formula  repre¬ 

sentation  using  "terrestrial"  gravity  anomaly  data  with  a  space  domain 
half-wavelength  resolution  of  0‘.5  55  km  on  the  surface  of  the  earth). 

Figure  5  presents  another  view  of  the  situation  by  plotting  from  the 
quantity  cml^r)  of  (3.31),  or  with  basically  the  same  numerical  values 

idso  from  cn^^n)  of  (3.35).  To  construct  Figure  5  a  given  altitude  was 
first  taken,  then  values  of  define  the  upper  limit  of  signal 

field  representations  were  searched  using  (3.31),  such  that  the  relative 
omission  error  CNmax  took  values  of  1%,  3X  and  lOX.  The  procedure  was 
repeated  for  altitudes  ranging  from  5  km  to  500  km.  For  example, 

Figure  5  says  that  for  an  altitude  of  30  km  we  need  a  representation  to 

at  least  n=(720,  540,  360)  in  order  to  keep  the  relative  omission  error 
down  to  (IX,  3X,  lOX).  Such  values  of  n  correspond  to  a  space  domain 
resolution  of  (15',  20',  30')  half-wavelength,  corresponding  to  about  (30, 
40,  55)  km  on  the  surface  of  the  earth.  Using  Figure  5  in  another  way, 
one  concludes  that  a  representation  to  n=360  incurs  less  than  IX 
relative  omission  error  for  altitudes  above  70  km,  incurs  3X  error  at  45 
km,  and  incurs  greater  than  lOX  error  for  altitudes  below  30  km. 

A  final  remark  should  be  made.  As  we  have  seen  the  radial  and 
horizontal  disturbance  signals  both  require  about  the  same  resolution  of 
representation  to  satisfy  certain  accuracy  levels.  For  example,  again 
from  Figure  5  we  see  that  at  an  altitude  of  5  km  and  desired  accuracy 
level  of  IX,  we  need  a  representation  to  about  n=2500.  This  corresponds 
to  a  need  in  an  integral  representation  for  5'x5'  gravity  anomalies 
(180*/2500^5' )  on  the  surface  of  the  earth.  This  is  true,  but  it  is  also 
true  that  the  5'x5'  values  need  be  given  only  within  a  very  limited  data 
cap  around  the  computation  point  and  not  all  over  the  earth.  This  is 
because  of  the  decay  of  the  integral  kernel  with  distance  from  the 
computation  point.  Farther  and  farther  away  from  the  computation 
point,  less  and  less  resolution  of  data  is  required.  In  this  aspect,  as 
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we  will  detail  in  the  next  section,  we  have  the  difference  that  the 
horizontal  disturbance  component  requires  significantly  larger  data  caps 
than  the  radial  component.  Therefore,  aside  from  the  spectral  study 
that  we  have  done  above,  we  also  need  a  space  domain  study  of  the 
actual  response  of  the  disturbance  signal  components  in  space  to 
terrestrial  gravity  anomaly  data  of  certain  resolution  and  distance  from 
the  computation  point.  The  spectral  study  shows  the  required  maximum 
resolutions,  but  the  space  domain  study  will  show  the  required  cap  sizes 
or  the  actual  spatial  location  and  resolution  of  the  required  Ag-data. 


3.3  Truncation  Theory 
3.3.1  Isotropic  Kernels 

Let  us  write  a  general  isotropic  integral  operation  on  the  gravity 
anomalies  Ag^i  producing  a  gravimetric  quantity  f,  as  follows: 


f  =  4;  J  I  K(t(-)  AgE  da  . 


(3.37) 


In  this  section  K(‘^)  is  the  isotropic  integral  operator  kernel,  that 
depends  on  the  angular  distance  f  but  not  on  the  azimuth  a  from  the 
computation  point  to  the  moving  integration  point.  The  coordinates  to 
use  in  (3.37)  are  the  geodetic  coordinates  (4,  X),  e.g.,  to  compute  the  i/ 
as  in  (3.17).  This  is  a  consequence  of  the  mapping  being  used  to 
convert  from  an  ellipsoid  to  a  spherical  boundary  surface  when  handling 
AgE-data  (Figure  2).  The  foregoing  discussions  can  easily  be  made  to 
apply  to  the  handling  of  boundary  values  Ag  that  refer  to  a  geocentric 
sphere  of  radius  rE  (those  defined  under  (3.9)),  by  just  changing  R  to 
rE>  ri  to  r,  4  to  ?,  and  AgE  to  Ag  in  the  formulas.  Examples  of  the  form 
(3.37)  are  the  expressions  (3.14)  for  T  and  (3.15)  for  We  can 

express  (3.37)  in  spectral  form  as  follows: 


f  =  2  YI.  ’ 

n=2 

consistent  with  the  expansions: 

K  =  YY.  Pn(y)  y  =  cos  i> 


4 

=  J  K(V)  P„(y)  dy 


(3.38) 


(3.39) 


(3.40) 


(3.41) 


4«e  =  l__  A«e 


n=2 


n 


AgEn  =  ^  J  J  AgE  Pn(y)  da  . 


(3.42) 


The  central  idea  in  truncation  theory  is  that  there  are  two  estimates 
that  we  denote  by  ^nd  Ag^*  available  for  Ag£,  these  estimates 

having  their  respective  errors  and  «*: 


AgE^  ~  Agf  +  (3.43) 

AgE*  =  Age  +  *•  .  (3.44) 


In  these  equations  the  superscript  "T”  denotes  the  estimate  from 
terrestrial  data  and  the  superscript  "s"  denotes  the  estimate  implied  by 
an  available  finite  set  of  spherical  harmonic  coefficients,  these  two  types 
of  estimates  being  the  ones  most  often  considered.  The  role  of 
truncation  theory  is  then  to  combine  the  information  from  the  above 
estimates  for  the  purpose  ^of  computing  the  unknown  quantity  f.  It  is 
generally  regarded  that  Agc^  will  be  used  as  carrier  of  high  frequency 
information,  while  Agf*  will  be  used  as  carrier  of  low  frequency 
information  about  the  gravity  field.  Considering  that  high  frequency 
information  is  associated  with  data  close  to  the  computation  point  while 
low  frequency  information  is  associated  with  data  farther  away,  then  the 
combination  of  Ag^^  and  Ag^*  to  produce  an  estimate  f  is  effected 
through  the  following  split  of  (3.37),  as  illustrated  in  Figure  6: 


f  =  f  +  Af 


(3.45) 


f  =  4^  J  J  (K-K»)  til  ^  I  *8tn 

n=2 


(3.46) 


Af  =  f  -  f  =  t;  J  f  (K-K»)  k-  e»  +  i  k»  AgE„ 


46 


where  we  have  the  kernel  for  the  spherical  harmonic  contribution  (see 
Figure  6): 


K* 


Rc  ,  0  ^  y,  4 

K  .  Vo  V  ^  " 


Vo  • • '  truncation  cap  radius 


(3.48) 


vrith  the  alternative  "modification”  functions  (Molodensky  et  al.,  1962; 
Meiaali  1971;  Wong  and  Gore,  1969;  Jekeli,  1980): 


K* 


Kx  =  0  [UnBodified  Molodensky  Truncation] 

K2  =  K(Vo)  [Meissl  Truncation] 
h 

kS  =  )  kn  Pn(y)  [Wong/Gore  Truncation] 

n=2 

h 

k5  =  )  S„  Pn(y)  ,  with  S„  such  that 

n=2 


(3.49) 


f  n  - 

J  (K-K4)^  sin  i>  -*  Bin.  [iBproved  Molodensky 
Vo  Truncation] 


In  (3.46)  and  (3.47)  the  kj;  are  the  coefficients  in  the  Legendre 
expansion  of  K*,  expressed  as  an  integral  of  K”  in  (3.40);  the  are  the 
surface  harmonics  of  the  error  c*,  expressed  as  an  integral  of  e*  as  in 
(3.42);  and  q  is  the  highest  degree  of  spherical  harmonic  representation. 
As  mentioned  under  (3.13)  of  the  last  section  we  do  not  use  in  the 

actual  computation  of  the  spherical  harmonic  coefficient  contribution  to  f 
(i.e.,  in  the  second  term  of  (3.46));  rather,  we  use  Ag^  as  given  by  (3.9) 
which  for  estimates  of  potential  coefficients  C*,,  becomes: 


r  (n-l) 


Y„„(?,  X)  . 


(3.50) 
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Kernel 


(truncation  cap  radius) 


Figure  6.  Splitting  of  the  Original  Kernel  K  for  the  Combination  of 
Terrestrial  Gravity  Anomaly  and  Spherical  Harmonic 


The  replacement  (3.50)  has  also  been  used  in  Rapp  (1981a,  second  term 
of  (34)). 


The  kernel  split  illustrated  in  Figure  6  conforms  to  the  idea  that 
detailed  (i.e.,  high  frequency)  terrestrial  information  need  be  used  only 
out  to  some  limited  radius  Vo  around  the  computation  point,  and 
therefore  the  kernel  (K-K")  for  the  terrestrial  contribution  is  set  to 
zero  outside  Vo>  From  the  space  domain  point  of  view  the  role  of  the 
modification  function  is  then  to  effect  a  certain  degree  of  continuity 
to  the  transition  from  the  non-zero  values  of  (K-K“)  for  V  <  Vo  to  its 
zero  values  for  V  ">  Vo>  equivalently,  the  modification  function 

provides  a  certain  degree  of  continuity  to  each  of  the  component 
kernels  (K-K*)  and  K*  of  the  total  kernel  K.  In  this  respect  the 
unmodified  Molodensky  truncation  (K‘==0)  does  not  provide  good 
continuity,  as  in  fact  in  this  case  the  kernels  (K-K*)  and  K*  become 
sharply  discontinuous  at  V  =  Vo*  Such  sharp  discontinuities  of  the 
kernels  in  the  space  domain,  both  in  function  values  and  derivatives, 
are  to  be  avoided  because  in  the  frequency  domain  they  directly  cause 
unwanted  oscillations  in  the  Legendre  coefficients  (k„  -  kj;)  and  k;;. 
These  oscillations  had  been  variously  called  "ringing”,  "ripples",  or 
"Gibbs  phenomenon"  in  Fourier  series  analysis.  If  we  have  these 
oscillations  then  we  are  not  getting  a  "clean"  separation  between  high 
and  low  frequency  information:  (K-K*)  may  have  significant  energy 
leaks  into  the  low  freqencies,  and  K*  may  have  significant  leaks  into  the 
high  frequencies.  In  terms  of  our  expression  (3.47)  for  the  total  error 
Af,  the  low  frequency  leaks  of  (K-K*)  will  enlarge  the  terrestrial  data 
propagation  error  (first  term  in  (3.47))  by  allowing  more  low  frequency 
energies  of  to  pass,  and  the  high  frequency  leaks  of  kj(  will  enlarge 
the  spherical  harmonic  coefficient  omission  error  (third  term  in  (3.47)) 
and  to  some  extent  the  spherical  harmonic  coefficient  commission  error 
(second  term  in  (3.47))  by  allowing  more  high  frequency  Agcn  ^nd  to 
pass. 


Several  authors  have  numerically  confirmed  the  significance  of 
carefully  choosing  the  modification  function  in  the  kernel  split  of  the 
type  shown  in  Figure  6.  We  mention  only  two  authors,  who  present 
several  comparisons  covering  all  the  famous  modification  functions  listed 
in  (3.49).  Jekeli  (1980)  presents  comparisons  for  the  case  of  the  Stokes 
kernel  and  Jekeli  (1982)  deals  with  the  Stokes,  Vening-Meinesz,  and 
Poisson  kernels.  Chen  (1982)  presents  comparisons  for  the 
Vening-Meinesz  kernel.  The  treatment  of  truncation  theory  for  the 
(non-isotropic)  Vening-Meinesz  kernel  is  only  slightly  different  from  an 
isotropic  case,  and  we  will  explain  this  below.  Also  below  we  will 
present  our  own  numerical  comparisons,  for  the  case  of  gravity 
disturbance  vector  computations  at  altitude.  In  these  computations  the 
basic  kernel  is  the  space-extended  Stokes  kernel,  which  is  differentiated 
to  yield  the  radial  disturbance  kernel  and  to  form  the  horizontal 
disturbance  kernel  (which  is  also  essentially  the  space-extended 
Vening-Meinesz  kernel). 
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The  decision  on  which  modification  function  is  the  best  is  usually 
based  on  least  RMS  value  of  the  total  error  (commission  plus  omission) 
of  the  spherical  harmonic  coefficient  contribution  to  f.  Specificallyi 
denoting  this  error  by  &f*,  we  have  from  (3.47): 


--it 


+  5  XI 

n=q+l 


(3.51) 


with  the  goal  to  have  the  least  value  for  the  globed  RMS,  expressed  as 
the  square  root  of: 


'(Af ) 


=  it 


(kS)* 


i  (k»)»  C„ 

n=q+l 


(3.52) 


Mere  the  ^C„  are  the  anomaly  error  degree  variances  implied  by  the 
spherical  harmonic  coefficient  errors  emd  the  C„  sure  the  anomaly  degree 
variances  themselves,  e.g.,  as  modeled  in  (3.36).  Equation  (3.52)  applies 
to  any  isotropic  kernel.  The  expression  for  the  non-isotropic  Vening- 
Meinesz  kernel  will  be  given  below.  To  summarize  the  numerical  compar¬ 
isons  (see  Jekeli,  1980,  1982;  and  Chen,  1982  for  details)  the  improved 
Molodensky  truncation  method  consistently  gives  better  results  (i.e., 
smaller  e(Af*))  than  any  of  the  other  truncation  methods  listed  in  (3.49). 
For  truncation  cap  radii  and  accuracy  levels  of  interest,  the  Meissl 
truncation  method  gives  better  results  than  the  unmodified  method  in 
the  case  of  the  Stokes  kernel,  but  yields  poorer  results  for  more  local 
kernels  (those  with  less  remote  zone  effects)  such  as  the  Poisson  and 
Vening-Meinesz  kernel.  The  Wong/Qore  method  is  only  of  interest  but 
has  not  been  specificcdly  designed  to  avoid  the  kernel  continuity  prob¬ 
lem  discussed  above,  and  in  fact  this  method  gives  significantly  poorer 
results  than  the  other  methods  (including  the  unmodified  method)  in 
typical  cases. 


3.3.2  Non-Isotropic  Kernels 


Let  us  now  turn  to  a  non-isotropic  integrsu  operation  on  the  gravity 
anomalies  Ag^.  Specifically,  we  will  discuss  the  integral  (3.16)  for  the 
computation  of  the  North-South  and  East-West  components  of  the  gravity 
disturbance  vector  in  space.  We  can  write  (3.16)  in  the  form: 


= d  a  ‘ft 


(3.53) 


I 

I 


where  Viii)  is  the  space-extended  Vening-Meinesz  function  as  given  in 
(3.20).  The  spectral  form  of  (3.53)  can  be  written  as  (see  (3.12)): 


(3.54) 


where  now  we  are  using  an  expansion  in  terms  of  the  associated 
Legendre  functions  of  the  first  order  (see  (3.28)): 


n=2 

.  cos?  *X 

V  =  ^  ^  P„,(y) 

n=2 

with  the  coefficients  (see,  e.g.,  Chen,  1982): 

1'  '(i*’  • 


(3.55) 


(3.56) 


Analogous  to  the  split  (3.45)-(3.47)  for  an  isotropic  integral  operation, 
we  now  have  the  following  split  of  (3.53)  for  the  computation  of 
estimates  3^  and  from  a  combination  of  terrestrial  gravity  anomaly 
and  spherical  harmonic  coefficient  information: 


(3.57) 


=  ^  J  J  (V-V-) 


cos  a 
sin  a 


tw 


cosV  *X 


(3.58) 


cos  a 
sin  a 


cos? 


IE' 

n=q+l 


cos?  *X 


(3.59) 


The  splitting  of  the  original  function  V  into  the  function  (V-V»)  for  the 
terrestrial  anomaly  contribution  and  the  function  V*  for  the  spherical 
harmonic  contribution  is  completely  analogous  to  the  splitting  of  K  into 
(K-K*)  and  K*  (see  Figure  6).  A  modification  function,  which  can  now 
be  denoted  as  is  also  used.  The  difference  is  that  now  all 

expansions  are  in  terms  of  P„t(7)  rather  than  ?„(/)»  e.g.,  in  such 
aquation  as  (3.49).  The  replacement  indicated  in  (3.50)  is  also  used. 
The  total  error  (commission  plus  omission)  of  spherical  harmonic 
coefficient  contribution  is  given  as  the  sum  of  the  second  and  third 
terms  in  (3.59): 


Adg 

f  ^  1 

4E«s 

» 

«? 

AdX 

n=2 

.  cos?  *X 

n=q+l 

cos?  •X 

(3.60) 


The  total  horizontal  error  from  (3.60)  is  obtained  from: 

(Atf„»)*  =  (4tff*)»  +  (A6x')*  • 


(3.61) 


The  goal  of  the  choice  of  modification  function  in  this  non-isotropic  case 
is  to  have  the  least  global  RMS  value  for  RMS  value  being 

expressed  as  the  square  root  of  (cf.  Heiskanen  and  Moritz,  1967,  pp. 
261-262): 


a*  (Ad 


n(n+l)(v«)*  dC„  + 


i  Y2  n(n+l)(vS)*  C„  . 

n=q+l 


no  •  •  •  .  n  .  •  o  •  o  .  -  .  •  .  •  -  ^  •  .. t a*'  o  ^  .  n  «  •  . n  o  ►  •  o  •  • 


(3.62) 


3»3.3  Truncation  Coefficients 

Now  let  us  discuss  how  to  obtain  the  so-called  truncation  coef¬ 
ficients  k|;  and  vg  needed  in  (3.46)  and  (3.58)  for  actual  estimations  and 
in  (3.42)  and  (3.62)  for  error  analyses.  We  will  discuss  the  coefficients 
for  the  computation  of  the  disturbance  vector  componentSi  for  which  the 
unsplit  kernels  are  (3.19)  for  the  radial  and  (3.20)  for  the  horizontal 
components.  The  starting  point  of  discussion  is  the  recursive  computa¬ 
tion  of  the  unmodified  truncation  coefficients  (i.e.,  those  corresponding 
to  the  use  of  K'==Ki‘^=0;  see  (3.49))  for  the  space-extended  Stokes  func¬ 
tion  given  in  (3.18).  Two  significantly  different  sets  of  recursive 
formulas  exist,  one  by  Shepperd  (1979)  and  the  other  by  Paul  (1983). 
Although  Paul's  recursives  are  claimed  to  be  more  stable  than 
Sheppard's  (see  Paul,  1983,  Table  3)  our  own  computations  showed  prac¬ 
tically  the  same  numerical  values  of  individual  truncation  coefficients 
from  the  two  solutions.  We  did  our  tests  (not  shown  here)  using  the 
program  given  in  Shepperd  (ibid.)  and  a  program  supplied  to  us  by 
Paul,  with  all  computations  implemented  in  double  precision.  We  used 
various  altitudes  and  compared  truncation  coefficients  from  degree  2  up 
to  some  reasonable  maximum  degree  of  required  field  representation  for 
the  altitude  in  question  (see  Figure  5).  Since  there  was  no  significant 
difference  between  using  Paul’s  or  Shepperd's  solution,  we  simply 
decided  to  use  Shepperd's  solution  for  all  subsequent  applications. 

Prom  the  unmodified  truncation  coefficients  for  the  space-extended 
Stokes  function,  Shepperd  (ibid.)  makes  extensions  to  compute  the  un¬ 
modified  truncation  coefficients  for  the  functions  (3.19)  and  (3.20)  of 
radial  and  horizontal  disturbance  computation,  and  gives  recursive 
formulas  and  programs  for  these  computations.  Shepperd’s  programs 
were  converted  to  double  precision  computations  for  our  applications. 

Now  given  the  unmodified  truncation  coefficients  for  the  functions 
(3.19)  and  (3.20)  of  radial  and  horizontal  disturbance  computation,  it  is 
relatively  easy  to  convert  these  coefficients  using  a  few  additional 
recursive  relations  to  obtain  the  truncation  coefficients  under  the 
Meissl,  Wong/Gore,  and  improved  Molodensky  truncation  methods.  For 
this  purpose  we  used  formulas  from  Jekeli  (1980),  which  were  given  for 
the  Stokes  function  but  are  actually  applicable  to  any  isotropic  kernel 
such  as  the  radial  disturbcmce  kernel  (3.19).  For  the  non-isotropic  case 
of  horizontal  disturbance  computation  we  used  formulas  from  Chen  (1982) 
for  the  conversion  from  the  unmodified  truncation  coefficients  from 
Shepperd’s  program  to  modified  truncation  coefficients  under  the  vari¬ 
ous  truncation  methods.  This  completes  the  discussion  of  how  to  obtain 
the  truncation  coefficients  for  radial  and  horizontal  disturbance 
computations. 


3.3.4  Numerical  Investigations 


For  error  analyses  using  (3.52)  and  (3.62)  we  needed  the  anomaly 
degree  variances  C„i  for  which  we  used  the  two-component  model  (3.36). 
We  also  needed  the  anomaly  error  degree  vcu'iances  AC„,  implied  by  the 
errors  in  the  spherical  harmonic  coefficient  set  to  degree  q.  Since  we 
wanted  to  test  several  reference  field  degrees  (q=20,  36,  180,  360)  we 
decided  to  use  a  model  for  6C„,  testing  two  models  A  and  B  derived 
from  a  white  noise  process: 


Model  A: 

6C(}  =  <ro  (2n+l)  ,  D  =  2,  3,  . . . ,  q 

with  <ro  such  that  SC*  =  0.5  Cq  .  (3.63) 

Model  B: 

tfC®  =0  ,  n  =  2,  3,  ....  10 

6C5  =  <fo(2n+l)  ,  n  =  11,  12,  . . . ,  q 

with  »o  such  that  6C*  =  0.3  Cq  .  (3.64) 

To  visualize  the  above  noise  processes  we  note  that  a  white  noise 
process  approjdmates  an  uncorrelated  error  process  in  blocks  of  certain 
size.  Specificedly,  if  we  have  uncorrelated  errors,  with  variance  m^(Sg), 
in  &g  in  blocks  of  size  9x6  (radians)  then  we  can  write  (Jekeli  and 
Rapp,  1980,  (10)): 


<ro(*n*(4«nin))  =  4;;  B*(4g). 


(3.65) 


From  (3.63)  and  (3.64)  we  can  also  express  as: 


(3.66) 


where  x  denotes  a  constant  factor  (x  =  0.5  in  (3.63),  and  x  =  0.3  in 
(3.64)).  A  harmonic  degree  q  corresponds  to  block  size  9  of: 


6  ■  -  (radians) 

q 


(3.67) 


Coobining  (3.65)  -  (3.67)  we  get: 


Usin^  (3.68)  and  the  two-component  model  (3.36)  for  Cq  we  can  now 
visualize  the  noise  process  in  (3.63)  (x=0.5)  as  beinf  uncorrelated  errors 
of  m(Ag)=(5,  6,  12,  14)  mgals  in  blocks  of  size  d=(9,  5,  1,  0.5)  degrees 
when  using  a  reference  field  to  degree  q=(20,  36,  180,  360).  For  (3.64) 
(x=0.3)  the  visualization  is  that  of  uncorrelated  errors  of  !n(Ag)=(3,  3.5, 
7,  8)  mgals  in  blocks  of  size  d=(9,  5,  1,  0.5)  degrees  when  using  a 
reference  field  to  degree  q=(20,  36,  180,  360).  In  (3.64)  there  is  an 
additional  removal  of  long-wavelength  errors  corresponding  to  harmonic 
degrees  2  to  10.  We  believe  that  (3.64)  is  a  more  realistic  portrayal  of 
available  data,  with  (3.63)  being  on  the  pessimistic  side. 

Since  only  the  improved  Molodensky  truncation  method  (with  tT=10, 
see  (3.49))  yielded  consistently  better  error  analysis  results  than  the 
unmodified  truncation  method,  we  do  not  present  results  for  the  Meissl 
and  Wong/Gore  truncation  methods  listed  in  (3.49).  Our  results  for  the 
improved  Molodensky  and  unmodified  Molodensky  truncation  methods  are 
given  as  Figures  7  to  14.  The  figures  show  the  RMS  truncation  error 
(in  mgals)  from  (3.52)  and  (3.62)  versus  the  truncation  cap  radius  Vo  (in 
degrees).  Several  reference  field  degrees  q,  denoted  as  nref  in  the 
figures  with  nref=20,  36,  180,  360,  were  tested  using  a  "high"  altitude  of 
100  km  and  a  "low"  edtitude  of  5  km.  Specifically,  the  Figures  7,  8,  9, 
10  give  results  for  the  radial  disturbance  case,  while  the  Figures  11,  12, 
13,  14  show  results  for  the  horizontal  case.  The  odd  numbered  Figures 
7,  9,  11,  13  contrast  the  case  of  having  commission  plus  omission  errors 
(solid  line)  against  the  case  of  having  omission  errors  only,  i.e.,  the 
case  of  perfect  reference  field  (circled  line),  under  the  unmodified 
Molodensky  truncation.  The  even  numbered  Figures  8,  10,  12,  14 
contrast  the  case  of  using  the  unmodified  Molodensky  truncation  (solid 
line)  agetinst  the  case  of  using  the  improved  Molodensky  truncation 
(circled  line).  The  upper  limit  n  of  the  summation  for  K'=  in  the 
improved  Molodensky  truncation  (see  (3.49))  was  taken  as  n=10,  which  is 
a  typical  value  that  can  be  chosen  (see  Jekeli,  1982).  A  value  of  h=20 
produced  slightly  larger  truncation  errors  for  a  given  cap  size.  Also, 
the  graphs  reflect  the  use  of  only  model  B  fur  (see  (3.64)).  The 

effect  of  using  model  A  for  6C„  (see  (3.63))  was  mainly  the  introduction 
of  an  additional  long  wavelength  error,  generally  making  the  graphs 
higher  by  about  0.5  mgal  for  small  cap  sizes  (’^o  <  3*)  decreasing  very 
slowly  to  0.2  mgal  for  large  cap  sizes  (Vo  up  to  40*  were  tested).  The 
magnitudes  of  RMS  truncation  errors  should  be  judged  against  the 
magnitudes  of  the  RMS  signals  themselves,  which  at  altitude  100  km  is 
about  20  mgals  and  at  altitude  5  km  is  about  36  mgals  for  both  the 
radial  and  total  horizontal  disturbances  (see  Figure  4). 

Looking  at  Figure  7  (solid  lines)  we  see  that  if  we  everytime 
tolerate  a  1%  truncation  error  (*0.2  mgal)  at  100  km  altitude,  then  we 
can  use  nref=20  (equivalent  to  the  use  of  mean  anomalies  in  block  sizes 
of  roughly  10*xl0*)  outside  a  cap  of  radius  V=13*  from  the  computation 
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point.  Inside  this  cap  Figure  7  shows  that  the  use  of  nref=20  will  make 
the  truncation  error  exceed  0.2  mgal,  and  in  order  to  avoid  this  we  can 
shift  to  the  use  of  nref=36  (^5*x5*  dg-data)  which  tfives  less  than  0.2 
mgal  error  down  to  For  ^  <  9*  we  can  shift  to  nref=180  (“1‘xl* 

dg-data)  and  this  will  control  the  error  to  below  0.2  mgal  down  to 
For  y=2*  to  3*  Figure  7  shows  that  nref=360  (^0!5x0t5  dg-data)  can  be 
used  to  keep  the  truncation  error  below  0.2  mgal.  For  ‘^  <  2*,  we  can 
extrapolate  from  Figure  7  that  in  order  to  keep  the  error  below  0.2  mgal 
we  can  use  higher  and  higher  order  reference  field  (denser  and  denser 
dg-data).  Alternatively,  for  ^  <  2*  we  can  use  a  more  accurate 

reference  field  to  degree  360  (more  accurate  OtSxO'.S  dg-data),  because  in 
fact  Figure  7  shows  that  an  errorless  reference  field  to  degree  360  is 
sufficient  to  define  the  radial  disturbance  at  altitude  100  km  with  a  1% 
relative  omission  error  (this  result  is  consistent  with  Figure  5  of  the 
last  section).  As  mentioned  at  the  end  of  the  last  section  and  now 
illustrated  by  Figure  7,  there  is  no  need  to  maintain  a  single  data 
density  all  over  the  earth.  Instead,  dg-data  of  less  and  less  resolution 
as  indicated  above  can  be  used  farther  and  farther  away  from  the 
computation  point.  Therefore,  graphs  from  truncation  theory  such  as 
those  we  show  here  directly  reveal  information  about  the  response  of 
gravimetric  quantities  to  terrestrial  gravity  anomaly  data  of  certain 
resolution  and  distance  from  the  computation  point. 

Comparing  Figure  11  with  Figure  7,  one  can  see  that  the  horizontal 
disturbance  requires  significantly  larger  data  caps  (more  than  2  times 
larger  than  the  radial  case)  in  order  to  keep  truncation  errors  below, 
e.g.,  a  0.2  mgal  level.  Therefore,  the  horizontal  disturbance  is  much 
more  sensitive  to  remote  zone  data  than  the  radial  disturbance.  This  is 
true  also  at  the  low  altitude  of  5  km,  where  Figure  13  (5  km)  shows 
only  slightly  different  curves  from  Figure  11  (100  km)  in  the  horizontal 
case,  while  Figure  9  (5  km)  shows  even  smaller  cap  requirements  than 
Figure  7  (100  km)  in  the  radial  case. 

To  summarize:  using  a  reasonably  typical  case,  assume  we  have 
detailed  gravity  anomaly  data  (e.g.,  5'x5')  in  a  cap  of  radius  3’  around 
the  computation  point.  In  quoting  truncation  theory  errors,  only  the 
errors  contributed  from  outside  the  inner  cap  are  considered,  with  the 
error  coming  from  inside  the  cap  being  propagated  separately  using  the 
first  terms  in  (3.47)  or  (3.59),  if  desired.  Now  in  order  to  account  for 
remote  zone  effects  (V  >  3*)  assume  that  we  use  a  reference  field  to 
degree  nref=180,  e.g.,  the  Rapp  (1981)  field.  Assume  further  that  the 
180-field  implies  the  errors  5C„®  given  by  (3.64)  (not  an  unreasonable 
assumption  for  the  Rapp- 180  field).  Then  the  RMS  truncation  error 
under  the  unmodified  Molodensky  truncation  theory  has  a  relative  vsdue 
of  under  IX  at  altitudes  100  km  and  5  km  for  the  radial  disturbance 
case,  and  can  reach  4X  at  100  km  and  2X  at  5  km  for  the  horizontal 
disturbance  case.  With  the  additional  presence  of  long  wavelength 
(n  *  10)  errors  in  the  180-field  and  implying  the  errors  given  by 

(3.63),  then  the  relative  truncation  errors  became  2X  at  altitude  100  km 
but  still  IX  at  5  km  in  the  radial  case,  and  may  reach  6X  for  altitudes 
100  km  and  5  km  in  the  horizontal  case. 


Looking  at  the  Figures  8,  10,  12,  14  contrasting  the  unmodified 
versus  the  improved  Molodensky  truncation  methods,  one  can  see  that 
the  improved  method  gives  better  results  for  relatively  large  cap  sizes 
only.  From  the  practical  point  of  view,  anticipating  the  use  of  nref=180, 
then  we  can  say  that  for  cap  sizes  and  accuracy  levels  of  interest  the 
improved  Molodensky  method  offers  no  significant  gain  over  the 
unmodified  case  in  both  radial  and  horizontal  disturbance  computations. 


4.  ACCOUNTING  FOR  THE  TOPOGRAPHY  IN  MODELING  THE  SPATIAL 
DISTURBANCE  VECTOR,  A  CONTINUOUS  APPROACH 

In  Chapter  1  we  introduced  the  gravity  anomaly  Ag,  related  to  the 
disturbing  potential  T  through  (1.9): 


4g 


1  i2 

y  *h 


T 


(4.1) 


or,  in  spherical  approxination,  through  (1.12): 


Ag  = 


(4.2) 


! 

If  the  given  gravity  anomalies  refer  to  an  equipotential  surface, 
then  the  modeling  of  the  anomalous  gravity  field  by  classical  integral 
formulas  (Heiskanen  and  Moritz,  1967,  sec.  6-4)  would  be  correct  to  the 
;  order  of  the  earth's  flattening  (ibid.,  pp.  240-241).  This  is  usually 

j  considered  a  sufficient  accuracy.  In  reality,  however  the  gravity 

!  anomaly  data  refer  to  the  earth’s  topography,  which  may  significantly 

deviate  from  being  an  equipotential  surface.  Therefore,  there  is  a 
problem  of  accounting  for  the  topography  in  the  modeling  procedures. 


In  this  chapter  we  will  be  concerned  with  the  above  topography 
problem.  We  will  present  a  composite  model  that  accounts  for  the 
topography  to  a  good  degree  of  approximation,  and  which  at  the  same 
time  is  convenient  for  practical  applications.  In  essence,  the  composite 
model  consists  of  the  spherical  harmonic  model  to  represent  the  long 
wavelength  part  of  the  field,  the  residual  topographic  model  to 
represent  the  short  wavelength  part  of  the  field,  and  the  classical 
integral  model  to  represent  the  remaining  (medium  wavelength)  part  of 
the  field.  In  the  next  chapter  we  will  additionally  examine  the  least 
squares  collocation  and  Bjerhammar’s  Dirac  type  models,  which  are  two 
other  models  that  can  be  introduced  to  account  for  the  topography  in 
modeling  the  external  gravity  field. 


4.1  Some  General  Considerations  on  the  Topography  Problem 


Some  general  characteristics  of  the  topography  problem  can  be 
concluded  by  considering  the  simple  crustal  density  model  used  also  by 
Moritz  (1968),  Lambeck  (1979),  Jekeli  (1981),  and  Rapp  (1982a).  Prom 
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Jekeli  (ibid.,  (2.130)}  we  have  that  the  spherical  harmonic  coefficients  hn„ 
of  the  equivalent  rock  topography  generate  the  spherical  harmonic 
coefficients  gj},,  of  what  we  will  call  "topographic  gravity  anomalies",  as 
follows: 


=  4.IKJP  (a^)  5^3  (4.3) 

where  R  is  a  mean  earth  radius  (6371  km),  G  is  the  constant  of 
gravitation,  and  p  is  a  constant  crustal  density  (2.67  g/cm’).  Equation 
(4.3)  assumes  that  the  topographic  masses  have  been  condensed  onto  the 
mean  earth  sphere  of  radius  R,  and  the  gj].  refer  to  this  same  sphere. 
Under  an  Airy-Heiskanen  isostatic  compensation  model,  compensating 
masses  at  depth  D  generate  spherical  harmonic  coefficients  g^M  what 
we  will  call  "isostatic  gravity  anomalies",  referred  to  the  mean  earth 
sphere  of  radius  R.  Assuming  that  the  compensating  masses  have  been 
condensed  onto  a  geocentric  sphere  of  radius  (R-D),  then  the  g^,,  are 
simply  the  negative  attenuated  version  of  gj}„,  as  follows: 


gSm  =  -(^)"  iHm  •  (4.4) 

The  coefficients  of  the  total  gravity  anomaly  referred  to  the  mean  earth 
sphere  can  be  expressed  as: 

inm  ~  Snm  gnm  •  (4.5) 

Combiniag  (4.3)'(4.5)  we  also  have: 


>'0 
.•v  ,'‘-' 


Given  the  anomaly  degree  variances  as  in  (3.36),  then  using  (4.6)  we 
can  propagate  €„  to  the  degree  variances  C|;  of  topographic  gravity 
anomalies,  as  follows: 


cS  = 


(4.7) 


Equation  (4.7)  expresses  the  degree  variances  of  the  gravity  anomaly 
field  generated  by  topographic  masses,  i.e.,  the  field  implied  by  the 
spherical  harmonic  coefficients  g!;„  of  (4.3).  The  presence  of  the 
isostatic  depth  D  in  the  equation  is  caused  by  the  following:  we  first  of 
all  used  D  to  define  the  (4.4),  then  with  the  condition 

(4.5)  that  the  sum  of  the  g{!„-  and  g^„-fielda  is  equal  to  the  observed 


field  g„„,  the  dependence  of  the  2{;„-field  on  the  depth  D  follows. 

Now  using  (4.4)  and  (4.7),  we  have  the  degree  variances  Cg  of  the 
gravity  anomalies  generated  by  isostatic  compensating  masses: 


fR-D] 

Cc  -  ^  «  *  C 

(‘  -  (¥1  ] 


(4.8) 


Using  the  degree  variances  C„,  C|;,  Cg  and  assuming  that  we  can 
computationally  isolate  the  gravitational  effect  of  topographic  masses  at 
harmonic  degrees  n  ^  N,  then  we  can  define  various  types  of  fields  with 
their  corresponding  gravity  anomaly  degree  variances,  as  shown  in 
Table  3. 

As  shown  in  Table  3  the  "total  field”  is  the  total  observed  field, 
with  associated  anomaly  degree  variances  C„  as  modeled  for  example  by 
(3.36)  for  n=2  to  -.  The  "residual  topographic  field”  is  generated  by 
the  gravitational  influence  of  topographic  masses  beyond  a  given 
harmonic  degree  N.  Specifically,  the  anomaly  degree  variances  for  this 
field  are  assumed  to  be  zero  for  n=2  to  N,  and  equal  to  the  Cf*  of  (4.7) 
for  n  >  N.  Analogously,  the  "residual  isoetatic  field"  is  generated  by 
the  gravitational  i'  *Tuence  of  isostatic  compensating  masses  beyond 
degree  N,  resulting  in  anomaly  degree  variances  of  zero  for  n=2  to  N 
and  the  of  (4.8)  for  n  >  N.  The  "residual-topography  reduced  field" 
is  found  by  subtracting  the  residual  topographic  field  from  the  total 
field,  leaving  a  field  with  anomaly  degree  variances  of  C„  for  n=2  to  N 
(as  in  the  total  field)  and  for  n  >  N  (as  in  the  residual  isostatic 
field).  Of  interest  also  is  the  "N-field”  (for  example,  called  the  180-field 
for  N=180)  which  is  the  part  of  the  total  field  that  is  implied  by 
harmonic  degrees  less  than  or  equal  to  N,  giving  degree  variances  of 
for  n=2  to  N  and  zero  for  n  >  N.  Note  that  the  N-field  theoretically 
results  by  subtracting  from  the  total  field  the  gravitational  influence  of 
both  the  residual  topographic  masses  and  their  compensation. 

For  first  numerical  considerations,  let  us  illustrate  that  the 
residual-topography  reduced  field  is  much  smoother  (i.e.,  with  much  less 
power  in  the  high  frequencies)  than  the  original  "total  field."  To  do 
this  for  the  disturbance  vector  field  at  altitude  the  construction  of 
Figures  4  and  5  of  the  previous  chapter  was  repeated,  this  time  using 
degree  variances  C„  for  n=2  to  180  and  for  n=181  to  5400,  instead  of 
using  Cp  entirely  for  n=2  to  5400  as  was  done  in  Figures  4  and  5.  For 
C„  the  two-component  model  (3.36)  was  again  used,  and  for  Cg  the 
standard  depth  of  compensation  D=30  km  was  assumed  in  (4.8).  For 
smoother  looking  graphs  around  n=180,  a  linear  transition  was 
arbitrarily  provided  between  C,*o  and  C^joo  thereby  avoiding  a  sharp 
discontinuity  that  originally  occurred  between  Ciao  and  Again, 

graphs  obtained  for  the  radial  and  horizontal  disturbance  components 
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Table  3.  Anomaly  Degree  Variances  of  Various  Types  of  Fields. 


Degree  Variances  on  a  Sphere  of 
Radius  R  =  6371  kn 


Field  Designation 


Total  Field 


Residual  Topographic 
Field 


Residual  Isostatic 
Field 


Residual-Topography 
Reduced  Field 


N-field  (e.g. , 
180-field  for  N=180) 


were  practically  identical  so  that  only  graphs  for  the  radial  component 
are  shown.  The  final  graphs  are  shown  as  Figures  15  and  16, 
corresponding  to  Figures  4  and  5,  respectively. 

For  example,  Figure  15  says  that  for  the  residual-topography 
reduced  field  the  use  of  a  radial  or  horizontal  disturbance 
representation  to  harmonic  degree  n=360  leaves  only  about  (1.8,  1.3,  0.7) 
mgals  to  be  resolved  at  (5,  10,  20)  km  altitude,  whereas  for  the  total 
field  (Figure  4)  these  values  are  correspondingly  (16.2,  10.2,  4.5)  mgals. 
On  the  other  hand.  Figure  16  shows  shows  for  example  that  for  the 
residual  topography-reduced  field,  at  an  altitude  of  30  km,  we  need  a 
representation  to  only  n=(390,  290,  195)  in  order  to  keep  the  relative 
omission  error  down  to  (IX,  3X,  lOX),  whereas  the  corresponding 
numbers  for  the  total  field  (Figure  5)  are  n=(720,  540,  360).  Even  at  the 
lowest  altitude  shown  (5  km)  the  reduced  field  needs  a  representation  to 
only  n=620  whereas  the  total  field  needs  up  to  n=2500  for  IX  relative 
omission  error.  Thus,  using  a  simple  crustal  model  we  see  in  a  general 
way  a  considerable  smoothing  of  the  field  by  removal  of  high-frequency 
topographic  effects;  real  data  computations  later  this  chapter  will 
further  illustrate  this  smoothing.  An  alternative  would  be  to  remove 
both  the  effects  of  high  frequency  topographic  masses  and  their 
compensation;  however,  we  did  not  choose  this  approach  because  of  some 
practical  reasons  stated  in  Section  4.2.3. 

A  second  numerical  investigation  of  interest  would  be  to  answer  the 
question:  how  sensitive  are  the  various  fields  of  Table  3  to  errors  in 
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the  vertical  location  of  the  boundary  surface  to  which  the  given 
anomalies  describing  the  field  refer?  The  vertical  location  is  defined  by 
the  heights  of  boundary  surface  points  above  or  below  the  ellipsoid  or 
the  geoid.  Numerical  values  for  such  sensitivities  are  important  because 
they  directly  indicate  the  errors  committed  in  the  classical  integration 
procedures  in  which  the  topography  of  the  earth  is  ignored;  in  such 
procedures  the  given  anomalies  that  actually  refer  to  the  physical 
surface  of  the  earth  are  arbitrarily  considered  to  refer  to  the  geoid  (or 
to  the  same  accuracy,  the  ellipsoid). 

In  considering  the  above  sensitivity  problem  let  C„  be  a  general 
symbol  to  denote  the  anomaly  degree  variances  (on  a  sphere  of  radius 
R)  of  any  of  the  fields  in  Table  3.  For  simplicity  we  will  assume  a 
constant  error  d  (e.g.,  d=1.5  km)  in  the  vertical  position  of  the 
boundary  surface  to  which  the  given  anomalies  refer.  Specifically,  we 
consider  that  the  correct  procedure  would  be  to  first  analytically 
downwcurd  continue  the  given  surface  anomalies  to  the  geoid  through  the 
distance  d,  and  then  to  use  the  downward  continued  values  in  a 
classical  integral  representation  of  the  external  gravity  field.  The 
erroneous  procedure  would  be  to  directly  use  the  surface  anomalies  in 
the  integral  representation  formulas  mistakenly  considering  that  these 
anomalies  refer  to  the  geoid.  The  anomaly  error  degree  variances  in 
this  case  would  be  the  degree  variances  of  the  difference  between  the 
downward  continued  anomalies  and  the  given  surface  anomalies.  We 
therefore  have  the  following  anomaly  error  degree  variances  for  a  field 
from  Table  3: 


in  which  the  second  term  inside  the  brackets  represents  the  appropriate 
factor  to  downward  continue  the  surface  anomalies  to  the  geoid,  through 
the  distance  d. 


We  should  note  that  there  are  upper  limits  to  the  values  of  harmonic 
degree  n  and  downward  continuation  depth  d  that  we  can  use  in  order 
for  (4.9)  to  keep  giving  numerically  meaningful  results.  Too  large 
values  of  either  or  both  n  and  d  would  cause  the  divergent  character  of 


harmonic  anedytical  downwEn-d  continuation  to  play  a  significant  role, 
yielding  meaningless  numerical  results.  Numerical  investigations  on  the 
divergence  of  downward  continuation  of  the  earth’s  gravity  field  can  be 
found,  for  example,  in  Jekeli  (1981).  Our  numerical  experience  with 
downward  continuation  (details  next  chapter)  is  in  the  context  of  the 
Bjerhammar  (1964,  1975)  approach  of  analytically  downward  continuing  a 
finite  set  of  surface  gravity  anomalies  down  to  an  internal  sphere. 
Based  on  getting  useful  results  while  downward  continuing  a  resolution 
of  5'x5'  values,  corresponding  to  down  to  depths  d  on  the 

order  of  5  km,  we  conclude  that  for  the  purposes  of  the  present 


67 


discussions  we  can  use  (4.9)  with  n„ax=2500  and  d=1.5  km,  as  we  will  do 
below. 

The  {global  RMS  (root  mean  square)  effect  of  on  computed  radial 
disturbances  at  altitude  can  be  obtained  by  using  AC„  in  place  of  Cp  in 
(3.29): 


(4.10) 


The  percent  error_  in  computed  radial  disturbance  due  to  the  error 
degree  variances  ACp  can  be  expressed  as: 


X  Error  =  x  lOOX  , 

<T{6r) 


(4.11) 


where  &<r(6r)  is  the  RMS  error  in  radial  disturbance,  obtained  as  the 
square  root  of  A<r*(6p)  in  (4.10);  and  <^(5^)  is  the  RMS  value  of  the  total 
radial  disturbance  signal  itself,  obtained  as  the  square  root  of  <r^(Ar)  in 
(3.29).  As  before,  there  is  no  need  to  discuss  the  case  of  the  horizontal 
disturbance,  as  practically  the  same  numerical  results  are  obtained  for 
both  the  radial  and  horizontal  cases. 

In  Table  4  we  show  the  RMS  errors  A<r((5r)  from  (4.10)  for  various 
types  of  fields  from  Table  3,  and  various  altitudes  of  disturbance 
computation.  For  each  altitude  Table  4  also  shows  the  total  RMS 
disturbance  signal  <r(Ar)  from  (3.29).  The  following  numerical  values 
were  used  to  construct  Table  4:  R  =  6371  km;  r,  =  R  +  altitude;  vertical 
position  error  d=1.5  km;  compensation  depth  D=30  km;  N=180; 
two-component  Cp  model  from  (3.36);  and  all  summations  to  infinity  where 
truncated  to  n=2500.  For  perhaps  easier  visualization  we  also  plotted 
the  X  error  from  (4.11)  as  Figure  17. 

It  is  seen  that  up  to  altitude  30  km  the  errors  caused  on  the 
residual  isostatic  field  and  the  180-field  may  still  be  considered 
significant,  but  are  considerably  smaller  than  the  errors  caused  on  both 
the  total  field  and  the  residual  topographic  field.  This  result  could 
have  been  expected  since  the  residual  isostatic  and  180-fields  are 
"smooth"  fields,  the  smoothness  being  directly  caused  by  the  isolation  of 
the  residual  topographic  field  (Figures  15  and  16).  For  smooth  fields 
the  downward  continued  anomalies  and  the  given  surface  anomalies  are 
relatively  close  in  magnitude,  leading  to  anomaly  degree  variances  ACp 
from  (4.9)  having  low  energy,  and  low  propagated  energy  through  (4.10), 
at  all  frequencies.  Aside  from  confirming  a  general  expectation,  Table  4 
and  Figure  17  are  useful  in  providing  a  feeling  for  actual  errors 
incurred  due  to  an  error  in  defining  the  vertical  position  of  the 
boundary  surface  to  which  the  anomaly  data  refer.  The  error 
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Table  4.  RMS  Error  in  Various  Types  of  Disturbance  Fields  Due  to  a  1.5 
km  Error  in  the  Vertical  Position  of  the  Boundary  Surface  for 
the  Anomaly  Data. 


Altitude 

5  km 

10  km 

30  km 

100  km 

500  km 

SMS  Disturbance 
of  Total  Field 
(mgal) 

36.92 

32.63 

26.04 

20.60 

13.11 

Field  Designation 

RMS 

Error  (mj 

;al) 

Residual  Topo¬ 
graphic  Field 

2.84 

1.75 

0.53 

0.04 

0.00 

Total  Field 

2.63 

1.57 

0.49 

0.11 

0.02 

Residual 

Isostatic  Field 

m 

0.35 

0.16 

0.02 

0.00 

180-field 

0.49 

0.44 

0.30 

0.10 

0.02 

0  1  2  3  i  ’  5  ’  6  r  8 
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Fi^re  17.  Percentage  Error  of  Various  Disturbance  Fields  at  a  Given 
Altitude,  Due  to  a  1.5  km  Error  in  the  Vertical  Position  of 
the  Boundary  Surface  for  the  Anomaly  Data. 


magnitudes  shown  are  in  close  agreement  with  real  data  computations,  to 
be  presented  later  this  chapter,  over  a  mountainous  area  in  New  Mexico. 


The  considerations  of  this  section  provide  the  rationale  for  a 
practical  set  of  procedures  to  be  detailed  in  the  next  sections,  for 
improving  the  classical  integral  representation  where  the  topography  of 
the  earth  is  neglected  in  modeling  the  external  gravity  field. 
Essentially,  detailed  (1  km  x  1  km)  topographic  height  information  and  a 
reference  topography  to  spherical  harmonic  degree  180  will  be  used, 
along  with  an  assumed  standard  crustal  density  (2.67  g/cm^),  to  isolate 
the  field  roughness  caused  by  the  residual  topographic  field  beyond 
degree  180.  A  spherical  harmonic  model  of  the  gravity  field  to  degree 
180  will  also  be  removed  to  minimize  remote  zone  effects.  The  remaining 
smooth  field,  which  is  primarily  the  residual  isostatic  field,  will  then  be 
represented  by  classical  integral  formulas.  The  smoothness  of  the  latter 
field  will  allow  the  minimization  of  the  errors  associated  with  the 
assumption  of  the  classical  integral  representation,  namely,  that  the 
boundary  surface  (which  in  reality  is  the  earth’s  topography)  is  a  level 
surface.  The  spherical  harmonic  gravity  model  will  be  separately 
evaluated  using  available  operational  programs  (Rapp,  1982b;  Rizos, 
1979).  The  residual  topographic  field  will  be  evaluated  by  direct 
integration  of  the  gravitational  effects  of  the  residual  topographic 
masses,  using  the  operational  program  of  Forsberg  (1984). 


4.2  Complementary  Models 

Here  we  summarize  the  spherical  harmonic,  residual  integral,  and 
residual  topographic  models  that  we  will  combine  for  the  computation  of 
the  gravity  disturbance  vector  in  space.  Pre-requisite  for  the 
combination  of  the  different  models  is  the  ability  to  compute  their 
implied  gravity  anomaly  values  on  the  earth’s  surface  where  the  originad 
anomaly  data  refer,  and  to  compute  their  implied  disturbance  components 
at  the  computation  points  in  space.  We  now  summarize  the  equations 
and  operational  programs  for  the  required  computations  using  the  three 
models. 


4.2.1  The  Spherical  Harmonic  Model 

On  the  earth's  surface,  the  implied  gravity  anomaly  is  conveniently 
written  in  the  truncated  Taylor  series  form: 


4g«(h,  ♦,  X)  =  Ag»(rE,  ▼,  X)  +  [  ^7^  (te.  ?,  X)  ]  h  ,  (4.12) 

where  (See  Figure  18): 

s  superscript  to  denote  the  spherical  harmonic  model 


h,  ♦,  X  geodetic  coordinates  of  the  surface  point 

rg,  V,  X  geocentric  coordinates  of  the  point  on  the  ellip¬ 

soid,  along  the  same  normal  as  the  surface  point 
(h,  ♦.  X). 


radial  derivative  of  the  spherical  harmonic  model. 


For  the  maximum  degree  of  spherical  harmonic  model  that  we  will  use  in 
this  report  (n„ax=^80)  the  higher  order  terms  omitted  in  (4.12)  can  often 
be  neglected.  This  was  found  from  numerical  testa  related  to  the  "Test 
3"  to  be  discussed  in  Section  4.3.4.  Higher  order  terms  can  be  used  if 
this  becomes  necessary.  Note  also  that  the  valid  approximation  is  made 
that  the  radial  derivative  can  be  used  instead  of  the  strictly 

required  derivative  <Ag*/<h  along  the  ellipsoidal  normal.  We  have: 


»  n  n+2 

A«»(rE,  ▼.  X)  =  ^  YZ  IZ  (77)”  WT,  X)  (4.13) 

n=2  m=-n 


•(re.  ^  IZ  ZI  ("■^2)(n-l)  X) 


where: 


n=2  m=-n 


kM  geocentric  gravitational  constant 

a  equatorial  radius 


(4.14) 


fully  normalized  potential  coefficients  referred  to 
those  of  the  reference  ellipsoid. 


The  advantage  of  the  form  (4.12)  is  that  the  quantities  (4.13)  and  (4.14) 
that  refer  to  the  point  (rg,  ?,  X)  can  be  evaluated  in  a  very  fast  way  on 
a  grid  of  constant  latitudes  and  longitudes.  For  these  «:valuations  we 
will  use  the  program  by  Rizos  (1979).  A  fast  evaluation  is  needed  since 
for  typical  applications  we  need  dense  values  of  4gB(h,  4,  X)  that  are  on 
a  regular  (4,  X)-grid  although  at  varying  heights  h.  For  example,  for 
the  tests  to  be  presented  later  this  chapter,  we  needed  a  5'x5'  grid  of 
values  covering  a  7*x9*  area  leading  to  84x108=9027  values;  such  a  large 
number  of  values  would  be  prohibitive  to  compute  if  we  were  to  evaluate 
the  spherical  harmonic  series  directly  at  the  individual  (h,  4,  X)  - 
positions. 


where  (r,  T,  \)  are  the  geocentric  coordinates  of  the  computation  point. 
Since  for  our  applications  we  needed  evaluations  only  at  isolated  points 
in  space  (not  on  a  regular  grid),  we  used  the  operational  program 
described  in  Rapp  (1982b)  to  evaluate  (4.15)-(4.17)  on  a  point  by  point 
basis.  For  fast  eveduations  of  (4.15)-(4.17)  on  a  regular  grid  the 
program  by  Rizos  (ibid.)  can  be  modified  (which  we  have  done,  and  the 
program  is  available). 


4.2.2  The  Residual  Integral  Model 

On  the  earth’s  surface,  the  implied  gravity  anomaly  of  this  model  is 
found  by  subtracting  from  the  originally  given  surface  anomaly 
4g(h,  4,  X)  the  anomalies  implied  by  other  complemetary  models.  With 
4g*(h,  4,  X)  being  the  spherical  harmonic  gravity  anomaly  from  (4.12), 
and  denoting  by  Ag^(h,  4,  X)  the  topographic  gravity  anomaly  to  be 
discussed  later,  then  we  have  the  surface  gravity  anomaly  of  the 
residtial  integral  model  as: 

Ago  =  4g  -  Ag»  -  Agt  ,  (4.18) 


where  the  superscript  D  denotes  the  classical  direct  integration  method 
(Heiskanen  and  Moritz,  1967,  Sec.  6-4)  where  AgO  will  be  used. 

At  altitude,  the  disturbance  components  of  the  residual  integral 
model  are  expressed  as  integrals  of  the  AgO  of  (4.18),  under  the 
assumption  that  the  AgO  refer  to  a  level  surface,  particularly  the  geoidal 
surface.  For  our  purposes,  we  consider  the  slight  modification  that  the 
AgO  might  be  considered  to  refer  to  a  level  surface  at  some  elevation  h,„ 
from  the  geoid  (See  Figure  18),  where  h,„  is  some  mean  elevation  of  the 
anomaly  data.  Note  that  in  these  discussions  there  is  no  need  to 
distinguish  between  the  geoid  and  the  ellipsoid  as  a  reference  surface 
for  heights.  The  effect  of  introducing  h^  is  to  change  the  vertical 
distance  between  the  computation  point  in  space  and  the  reference  level 
for  the  data,  from  H  (when  h„=0)  to  Ho=H-h„  (when  h„^0),  where  H  is 
the  altitude  of  the  computation  point  above  the  ellipsoid  (or  the  geoid)  - 
see  Figure  18.  It  is  important  to  distinguish  between  Ho  and  H,  as  their 
relative  difference  can  be  large,  for  example,  H=5  km  and  h„=1.5  km  give 
Ho=3.S  km.  On  the  other  hand,  there  is  no  need  to  modify  the  mean 
earth  radius  R=6371  km  to  R'=R+h„,  as  the  relative  difference  between  R 
and  R'  is  small  (e.g.,  even  with  h,„=4  km,  h„/R  is  only  6  x  lO"*). 

We  now  summarize  the  expressions  for  the  disturbance  components 
at  altitude,  based  on  the  residual  integral  model.  For  the  implementation 
of  these  expressions  the  operational  program  fully  described  in  Rapp 
(1966b)  can  be  used.  From  Rapp  (ibid.)  the  disturbance  components  at 
geodetic  coordinates  (H,  4,  X)  become: 
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tff  =  ^  J  J  Ag®  Fa  cos  a  d<r 


fix  =  ^  J  J  Ag®  Fa  sin  a  d<r 


where: 


I  t»  [  ^  +  J  +  l-6D-t  cosf  (l3  +  6  In 


»  ^t..ri.3  *  3  fD+l-t  cosi( 

r.  =  f  simt  [  55  .  5  -  4  »  J  I  0 


D+l-t  cos 


H  +  Ho 


R  =  6371  km 


Ho  =  H  -  h„  (Figure  18) 


D*  =  (1  -  2t  cos  V  +  t*) 


(4.24) 


(4.25) 


(4.26) 


dff  =  cos  t  dt  dX  (rectangular  blocks  as  a  97^ 

integration  elements)  ^  *  ' 

angular  distance  between  the  projected  computation  point 
(R,  tp,  Xp)  and  the  moving  integration  point  (R,  X(,|) 

a:  azimuth  (clockwise  from  North)  from  the  projected  computation 

point  (R,  tp,  Xp)  to  the  moving  integration  point  (R,  ^|if,  X^). 

For  our  tests  we  used  a  modified  program  extracted  from  the  program  in 
Rapp  (ibid.),  where  the  modification  consists  of  a  simplified  input/output 
and  the  use  of  only  one  anomaly  block  size  (5'x5')  instead  of  all  the 
mean  anomady  block  sizes  2'.5  x  2'.5,  5'x5',  30'x30',  and  l*xl*  used  by 
the  original  program.  The  subtraction  of  Ag*^  in  (4.18)  removed  much  of 
the  roughness  of  the  Ag-field,  and  the  subtraction  of  Ag"  drastically 


m 
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reduced  remote  zone  effects,  so  that  for  the  purposes  of  our  tests  we 
used  only  5'x5'  mean  blocks  covering  7*x9*  area  around  the  computation 
points  in  New  Mexico.  To  be  able  to  make  computations  at  low  altitudes 
(values  of  Ho  down  to  3.5  km  were  used)  the  integrated  kernel 
evaluation  implemented  in  Rapp  (ibid.)  for  2'.5x2r.5  blocks  was  also 
implemented  in  the  modified  program,  at  least  for  those  5'x5'  blocks 
whose  centers  were  within  a  distance  of  10'  from  a  computation  point. 


4.2.3  The  Residual  Topographic  Model 

This  model  is  advocated  in  Porsberg  and  Tscherning  (1981)  and 
Forsberg  (1984)  for  use  in  general  gravity  field  modeling  problems.  The 
model  is  generated  as  the  gravitational  influence  of  residual  topographic 
masses  that  lie  between  the  actual  topography  with  elevation  h  and  a 
reference  topography  with  elevation  h^.  The  residual  masses  may  be 
positive  (h>h*)  or  negative  (h®>h).  The  reference  topography  is 
conveniently  a  spherical  harmonic  expansion  of  the  topography,  e.g.,  to 
degree  180,  but  it  can  also  be  any  averaged  version  of  detailed 
topographic  information. 

An  alternative  model  is  the  topographic/isostatic  model  used  in 
Sunkel  (1983a)  and  Moritz  (1983).  This  model  is  generated  as  the 
gravitational  influence  of  the  topographic  masses  (masses  between  the 
geoid  and  the  topography)  and  its  isostatic  compensation.  However,  for 
our  applications  we  choose  the  residual  topographic  model  because  of  its 
definite  computational  advantages  (Forsberg,  1984)  over  th 
topographic/isostatic  model,  namely,  less  remote  zone  and  edge  effects, 
no  necessity  to  compute  for  the  effects  of  compensating  masses,  and 
availability  of  a  simple  approximate  expression  for  the  model  gravity 
anomalies  on  the  earth's  surface  (see  (4.28)  below).  The  use  of  the 
residual  topographic  model  compares  favorably  with  the  use  of  the 
topographic/isostatic  model  as  complementary  model  to  collocation  for  the 
prediction  of  gravity  anomalies  and  deflections  of  the  vertical  in  a 
mountainous  test  area  in  New  Mexico  (Forsberg  and  Tscherning,  ibid.). 

Details  of  the  generation  of  various  gravimetric  quantities  (radial 
gravity  disturbance,  deflections  of  the  vertical,  height  anomaly,  and 
gravity  anomaly)  implied  by  the  residual  topographic  model  are  given 
along  with  an  operational  computer  program  in  Forsberg  (1984).  The 
quantities  can  be  computed  on  the  earth’s  surface  and  in  space  by 
integrating  the  gravitational  effects  of  the  residual  masses  of  assumed 
constant  density  (2.67  g/cm’  is  a  standard  crustal  density).  The  prism 
is  the  basic  integration  element  used.  Close  to  the  computation  point 
exact  prism  formulas  are  used,  while  farther  away  approximate  formulas 
from  spherical  harmonics  and  point  mass  effects  are  used  for  economy 
and  stability  of  computations. 


On  the  earth’s  surface,  the  implied  gravity  anomaly  4g^  of  the 
residual  topographic  model  can  be  obtained  either  by  direct  prism 
integration  using  Forsberg’s  program,  or  by  using  the  approximate 
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relation  (Forsberg,  ibid.t  p.  39): 


Agt  =  2»iGp  (h-h«)  -  tc  , 


(4.28) 


where: 


gravitational  constant 
density 

elevation  of  surface  point 
elevation  of  reference  topography 


terrain  correction. 


Equation  (4.28)  is  useful  in  the  case  when  the  tc  is  already  given  from 
the  data  records,  or  as  we  will  see  later,  in  the  case  when  tc  can  be 
neglected  due  to  Ag*^  being  sufficiently  far  from  the  spatial  points  where 
the  disturbance  components  are  being  computed.  In  the  latter  case  the 
approximate  residual  topographic  anomaly  is  given  by: 


Ag‘  =  27tGp  (h-h«) 


(4.29) 


At  altitude,  the  disturbance  components  implied  by  the 

residual  topographic  model  can  be  obtained  by  direct  prism  integration 
using  Forsberg’s  program.  A  slight  modification  is  necessary,  to 
compute  the  meridional  and  prime  vertical  disturbance 

components  instead  of  the  deflections  of  the  vertical;  we  have  (Rapp, 
1982b,  p.  6): 


=  -r  f 


fix  =  ~7  n  . 


(4.30) 


where  (  is  the  meridional  and  ti  is  the  prime  vertical  component  of  the 
deflection  of  the  vertical,  and  y  is  normal  gravity  on  the  ellipsoid 
linearly  attenuated  (0.3086  mgal/m)  to  account  for  the  height  of 
computation  point  above  the  ellipsoid  (Rapp,  ibid.,  p.  8). 


4.3  Numerical  Inveatiffaiiona 
4.3.1  The  Data 

For  our  numerical  tests  we  used  the  5"x5'  mean  values  of  surface 
gravity  anomaly  and  elevation  developed  over  a  7*x9*  area  in  New 
Mexico  in  Cruz  and  Laskowski  (1984).  For  certain  inner  zone 
computations  required  in  Forsberg’s  (1984)  prism  integration  program  we 
also  used  the  detailed  1  km  x  1  km  grid  point  elevation  data  supplied 
earlier  (April,  1983)  by  the  National  Geodetic  Survey  (NGS).  The 
reference  topography  used  in  the  residual  topographic  model  was 
generated  from  a  set  of  spherical  harmonic  coefficients  to  degree  180 
(tape  GS140,  file  15).  The  spherical  harmonic  gravity  model  used  was 
that  developed  in  Rapp  (1981b)  up  to  degree  180,  which  here  will  be 
referred  to  as  the  Rapp-180  field. 

For  immediate  visualization  Figure  19  and  20  show,  respectively,  a 
contour  map  of  5'x5'  mean  elevations  (contour  interval  =  50  meters)  and 
5'x5'  mean  gravity  anomalies  (contour  interval  =  5  mgfids)  in  a  portion 
of  our  New  Mexico  test  site.  The  anomalies  shown  are  actually  Faye 
anomalies,  i.e.,  terrain-corrected  anomalies  (ibid..  Figure  14).  The  Faye 
anomaly  is  expectedly  even  more  correlated  with  elevation  than  the 
terrain-uncorrected  free-air  anooaly  (Moritz,  1968,  p.  33).  This 
correlation  is  clearly  observed  from  Figures  19  and  20,  where  the 
contour  intervals  were  chosen  to  be  in  a  1  m  to  0.1  mgal  ratio  to 
approximate  the  standard  correlation  factor  of  0.1119  mgal/m. 

If  now  we  remcve  the  influence  of  tox>ographic  masses  by 
subtracting  the  Ag^  of  (4.28)  from  the  free-air  anomaly  or,  the  same 
thing,  by  subtracting  the  first  term  of  (4.28)  from  the  Faye  anomalies, 
then  we  should  see  a  considerable  smoothing  of  the  anomaly  field.  This 
is  shown  in  Figures  21  and  22.  In  Figure  21  the  entire  topographic 
masses  (i.e.,  h"=0  in  (4.28))  were  removed,  while  in  Figure  22  only  the 
residual  masses  (h”=reference  topography  to  degree  180)  were  removed. 
The  field  shown  in  Figure  21,  associated  with  the  subtraction  of  (2rrGph 
-  tc)  from  the  free-air  anomalies,  is  called  refined  Bouguer  anomalies  in 
Heiskanen  and  Moritz  (1967,  Sec.  3-3).  Figure  22  additionally  reflects 
the  removal  of  the  Rapp-180  field  as  Ag*,  so  that  the  figure  actually 
shows  the  residual  anomaly  Ag°  of  (4.18).  In  effect.  Figure  22  can  be 
obtained  by  subtracting  from  Figure  21  the  quantity  (Ag*  -  2rrGph*). 
Since  the  latter  quantity  can  be  interpreted  as  a  reference  Bouguer 
anomaly  (see  ibid.),  then  we  can  call  the  field  shown  in  Figure  22  as  a 
residual  refined  Bouguer  anomaly  field. 

Figure  21  is  indeed  much  smoother  than  Figure  20,  although  one  can 
observe  the  well-known  large  negative  bias  of  Bouguer  anomalies 
associated  with  the  removal  of  the  entire  topography  but  not  its 
isostatic  compensation.  Figure  22  avoids  the  bias  because  positive  and 
negative  masses  were  removed  or,  as  an  alternative  interpretation, 
because  relatively  high  frequency  topographic  masses  (beyond  degree 
180)  which  have  only  small  isostatic  compensation  were  removed. 
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Since  Figure  22  reflects  the  removal  of  the  lower  180  harmonics  and 
the  influence  of  topographic  masses  beyond  degree  180,  the  field  shown 
can  be  interpreted  as  the  residual  isostatic  field  discussed  in  Section 
4.1  (See  Table  4.1).  In  reality,  Figure  22  also  contains  the  errors  of  the 
4g*  and  4g^  fields  (see  (4.18)).  Nevertheless,  one  can  see  the  degree  of 
smoothness  of  Figure  22  compared  with  Figure  20.  In  fact  the  RMS 
anomaly  value  reduced  by  42%,  from  26  mgal  for  the  field  shown  in 
Figure  20  to  15  mgal  for  the  field  shown  in  Figure  22,  considering  all 
anomaly  values  in  the  7*x9*  test  area  described  in  Section  4.3.2.  This 
provides  a  good  real  data  confirmation  of  our  earlier  general  study  that 
the  topography-unreduced  field  (Figure  20)  might  require  a 
representation  to  degree  n=2S00  (5'x5'  resolution)  for  low-altitude 

disturbance  computations,  whereas  the  topography  reduced  field  might 
need  only  a  representation  to  n=620  (17'xl7'  resolution)  -  see  the 
discussions  related  to  Figures  15  and  16. 


4.3.2  The  Teat  Area 

For  testing  the  combination  of  the  various  models  discussed  in 
Section  4.2,  we  used  gravity  and  elevation  data  in  the  7*x9*  Now  Mexico 
area  shown  in  Figure  23.  The  figure  shows  the  reference  topography  to 
degree  180  with  a  contour  interval  of  50  meters.  The  figure  also  shows 
five  symmetrically  arranged  points  in  a  l*xl*  area,  these  points  being 
foot  points  of  five  vertical  trajectories  (i.e.,  lines)  along  which  test 
computations  were  performed.  Test  computations  were  performed  at 
points  along  the  trajectories,  at  altitudes  5,  10,  30,  100,  and  500  km 
above  the  ellipsoid.  The  locations  of  the  five  trajectories  were  chosen 
in  an  area  of  significant  gradient  of  the  disturbance  component  fields. 
To  show  the  gradients,  the  3*x3*  area  marked  off  in  Figure  23  was 
taken  for  plotting  the  disturbance  components  at  altitudes  0  km  and  30 
km  using  a  contour  interval  of  5  mgals.  The  plots  are  shown  in  Figure 
24.  The  disturbance  components  were  plotted  from  Ot25xO’.25  grid  point 
values  generated  from  the  Rapp-180  field  using  the  modified  Rizos 
program  mentioned  in  Section  4.2. 

Figure  25  shows  the  labeling  that  we  will  use  to  identify  the  test 
points  at  a  given  altitude.  Also  shown  is  the  area  (Area  I)  which 
contains  the  5'x5'  blocks  where  the  integrated  kernel  form  (Rapp, 
1966b)  was  used  in  the  numerical  integration  of  the  integral  model.  The 
integrated  kernel  form  avoids  the  large  errors  of  the  center-point 
kernel  evaluation,  for  integration  elements  close  to  the  computation  point 
at  low  altitudes.  Figure  25  also  shows  the  area  (Area  II)  where  rigorous 
values  of  Ag*  from  Forsberg’s  program  were  used  in  the  residual 
topographic  model.  Outside  Area  II,  it  will  be  shown  in  the  tests  that  it 
is  sufficient  to  approximate  by  the  first  term  of  (4.28). 
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4.3.3  Description  of  Teats  and  Numerical  Results 


Below  we  describe  the  various  tests  conducted  using  the  five 
vertical  trajectories  through  the  points  A,  Bt  C,  D,  E  shown  in  Figure 
25.  Again,  the  altitudes  tested  were  5,  10,  30,  100,  and  500  km.  To 
simplify  the  numerical  presentation  we  give  for  each  particular  test  and 
particular  disturbance  component,  the  numerical  results  for  only  one 
chosen  trajectory  omitting  the  results  for  the  remaining  four 
trajectories.  The  one  trajectory  to  be  presented  will  be  the  one 
showing  the  greatest  sensitivity  to  the  particular  test,  that  is,  the  one 
showing  the  largest  absolute  values  for  numerical  results.  We  now 
describe  the  various  tests,  with  the  results  for  each  presentation 
trajectory  being  shown  in  Table  5  at  the  end  of  the  test  descriptions. 

Test  1.  Remote  Zone  Error  of  Direct  Integration  Method  for  n=2  to  180 

The  Rapp-180  field  was  first  evaluated  to  provide  "true"  values  of 
disturbance  components  in  space  using  (4.15)-(4.17),  and  to  provide 
5'x5'  "data"  values  of  anomalies  on  the  ellipsoid,  in  the  7*x9*  test 
area,  using  (4.13).  The  data  values  were  then  input  into  the  direct 
integration  method  using  (4.19)-(4.21)  to  yield  computed  values  in 
space.  The  true  values  were  subtracted  from  the  computed  values 
to  yield  the  errors  for  this  test.  The  differences  are  caused  by  the 
neglect  of  the  anomaly  data  outside  the  data  area. 


Teat  2.  Remote  Zone  Error  of  Direct  Integration  Method  for  n=21  to  180 

Test  1  was  repeated  this  time  using  only  harmonics  from  n=21  to 
180. 


Teat  3.  First  Derivative  Effect  for  the  180-Pield 

5'x5'  values  of  the  second  term  in  (4.12)  were  first  computed  in  the 
7*x9*  test  area  using  the  Rapp-180  field  and  the  5'x5'  mean  values 
of  elevations  (not  the  reference  elevations  to  180).  The  resulting 
values  were  then  input  into  the  direct  integration  method  using 
(4.19)-(4.21)  to  yield  the  first  derivative  effects  for  this  test. 


Test  4.  Sensitivity  of  Field  Beyond  Degree  180  to  Height  Error 

The  spherical  harmonic  gravity  anomaly  4g*  from  the  Rapp-180  field 
and  (4.12)  was  first  subtracted  from  the  originally  given  surface 
anomalies  Ag,  to  form  5'x5'  values  of  (Ag-Ag“)  in  the  7*x9*  test 
area.  These  values  were  then  input  into  the  direct  integration 
method  using  (4.19)-(4.21),  first  with  h„=0  then  with  hm=1.5  km  in 
(4.25),  where  h„=1.5  km  was  an  approximate  mean  elevation  in  the 
test  area  (see  Figure  23).  The  sensitivities  for  this  test  were  then 
formed  by  subtracting  the  disturbance  results  for  the  case  when 
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h„=0  from  those  for  the  case  when  h„=l.S  km. 

Test  5.  Residual  Topographic  Model  Outer  Zone  Effects 

5'x5'  mean  elevations  and  30"x30"  grid  point  elevations  of  the  actual 
topograph^i  along  with  the  reference  topography  to  180f  were  input 
into  Forsberg’s  prism  integration  program  to  compute  the 
disturbance  components  4^^,  at  altitude.  The  standard 

crustal  density  of  2.67  g/cm’  was  used.  Two  versions  were 
computed,  one  using  data  in  the  entire  7*x9*  test  area,  and  the 
other  using  only  data  inside  a  2*x2*  area  centered  at  the 
computation  points  (Area  II  in  Figure  25).  The  results  for  the  2*x2* 
area  were  subtracted  from  the  results  for  the  7*x9*  area  to  yield 
the  desired  outer  zone  effects  (disturbances)  for  this  test. 

Test  6.  Error  of  Using  Ag*  =  2TtGg(h-h»)  as  a  Simplified  Residual 
Topographic  Anoamlv  in  the  Outer  Zone 

S'xS'  mean  values  of  Ag^=2TrGp(h-h*)  were  first  generated  in  the 
outer  zone,  outside  the  Area  II  of  Figure  25  but  inside  the  7*x9* 
test  area,  using  5'x5'  values  of  elevations  h  and  reference 
elevations  h*  to  degree  180,  and  using  2tiGpx0.1119  mgsl/m  associated 
with  p=2.67  g/cm’*  These  data  in  the  outer  zone  were  then  input 
into  the  direct  integration  method  (4.19)-(4.21),  with  h^^LS  km  in 
(4.25),  to  yield  the  spatial  disturbances  implied  by  the  outer  zone 
Ag^-data.  From  these  disturbances,  the  rigorous  outer  zone 
disturbances  from  Test  5  were  subtracted  to  yield  the  errors  for 
this  test. 


Test  7.  Sensitivity  of  the  Residual  Topographic  Field  to  Height  Error 

Rigorous  values  of  5'x5'  residual  topographic  anomalies  Ag^  on  the 
earth's  surface  were  first  generated  inside  the  2*x2*  inner  zone. 
Area  II  in  Figure  25,  by  direct  prism  integration  of  residual 
topographic  masses  using  Forsberg’s  program  (see  comment  above 
(4.28)).  Again,  the  input  data  were  5'x5'  mean  elevations  and 
30"x30"  grid  point  elevations  of  the  actual  topography,  reference 
topography  to  degree  180,  and  p=2.67  g/cm’.  These  Ag^  values  in 
the  inner  zone  were  combined  with  the  Ag^  values  developed  for  the 
outer  zone  in  Test  6.  The  combined  5'x5'  data  in  the  7*x9*  area 
were  then  input  into  the  direct  integration  method  (4.19)-(4.21),  first 
with  h„=0  then  with  h„=1.5  km  in  (4.25),  where  h„=1.5  km  was  an 
approximate  mean  elevation  in  the  test  area.  The  sensitivities  for 
this  test  were  then  formed  by  subtracting  the  disturbance  results 
for  h„=0  from  those  for  hm=1.5  km. 
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Teat  8.  Sensitivity  of  the  Residual  Integral  Field  to  Height  Error 


5'x5'  values  of  the  residual  integral  model  anomalies  Ag*’  from  (4.18) 
were  first  formed  in  the  7*x9*  test  area,  by  subtracting  from  the 
original  surface  anomalies  Ag  the  spherical  harmonic  anomalies  &g* 
from  (4.12)  and  the  residual  topographic  anomalies  &g^  (inner  zone) 
and  Ag^  (outer  zone)  from  Teat  7.  Again,  for  computing  Ag”  the 
Rapp-180  field  and  5'x5'  mean  elevations  of  the  actual  topography 
were  used.  The  Ag^^-data  were  then  input  into  the  direct  inte¬ 
gration  method  (4.19)-(4.21),  first  with  h„,=0  then  with  h^i^liS  km  in 
(4.25),  where  h,p=1.5  km  was  an  approximate  mean  elevation  in  the 
test  area.  The  sensitivities  for  this  test  were  then  formed  by 
subtracting  the  disturbance  results  for  h„=0  from  those  for  h„=1.5 
km. 


Teat  9.  Error  of  the  Direct  Integration  Method  in  Renresentint 
Residual  Topographic  Field 


5'xS'  values  of  residual  topographic  anomalies  on  the  earth's 
surface  were  first  formed  in  the  7*x9*  test  area,  from  the  Ag^ 
(inner  zone)  and  Ag^  (outer  zone)  values  from  Test  7.  These  values 
were  then  input  into  the  direct  integration  method  (4.19)-(4.21),  with 
hm=1.5  km  in  (4.25),  to  yield  the  spatial  disturbances 
(ix^)°  implied  by  the  integral  representation  method.  Corresponding 
rigorous  vedues  of  the  disturbances  were  generated 

from  Porsberg's  prism  integration  program,  using  5'x5'  mean 
elevations  and  30"x30"  grid  point  elevations  of  the  actual 
topography  along  with  reference  topography  to  degree  180  and 
p=2.67  g/cm’.  The  errors  for  this  test  were  then  as  the  differences 
[(A‘)0  -  an.  [(«4M°  -  AaM,  and  ((fixM“  -  ^xM- 


4.3.4  Discussion  of  Numerical  Results 

Test  1  shows  that  the  remote  zone  errors  caused  by  the  omission  of 
anomaly  data  outside  the  7*x9*  test  area  are  significant. 

Comparison  of  Tests  1  and  2  shows  that  a  large  part  of  the  remote 
zone  errors  are  caused  by  very  long  wavelength  harmonics  between 
degrees  2  to  20,  with  the  errors  caused  by  degrees  21  to  180  being 
below  1  mgal  but  can  still  be  considered  to  be  significant  biases.  The 
characterization  of  the  errors  as  biases  comes  from  the  fact  that  the 
errors  attenuate  very  slowly  with  altitude.  Extrapolating  the  results  of 
Tests  1  and  2,  it  can  be  expected  that  the  remote  zone  errors  caused  by 
harmonics  above  degree  180  are  negligible.  This  expectation  also  follows 
from  the  analysis  of  truncation  errors  as  presented  in  Chapter  3.  This 
means  that  it  will  be  sufficient  to  use  a  spherical  harmonic  model  to 
degree  180  (as  we  did  in  our  tests)  as  a  complementary  model  to  account 
for  remote  zone  data  outside  the  7*x9*  test  area. 
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Table  5.  Maximal  Values  (mgala)  of  the  Different  Quantities  Being 
Computed  in  Tests  1  to  9  of  Chapter  4  (see  text). 


30  ka 


100  ka 


500  ka 


-1.65 


-0.20 


-1.96 

-2.98 


Test  3  shows  that  for  the  purpose  of  subtracting  the  spherical 
harmonic  model  to  degree  180  from  surface  gravity  anomaly  data,  it  is 
not  sufficient  to  simply  use  the  values  of  the  spherical  harmonic  model 
on  the  surface  of  the  ellipsoid.  Rather,  the  first  radial  derivative  of 
the  field  has  to  be  used  also  as  in  (4.12)  to  account  for  the  varying 
heights  of  the  surface  points  above  the  ellipsoid.  Tests  not  shown  here 
indicate  that  the  effects  of  the  second  radial  derivative  are  one  order  of 
magnitude  below  those  of  the  first  radial  derivative,  and  thus  such 
effects  can  be  neglected.  Again,  we  point  out  that  the  evaluation  of  the 
spherical  harmonic  model  using  the  truncated  Taylor  series  (4.12)  is  a 
means  to  control  the  amount  of  computer  time  required  for  evaluations 
of  the  model  on  a  dense  grid  of  surface  i>oints. 

Test  4  shows  that  after  subtracting  the  spherical  harmonic  model  to 
degree  180  from  the  original  field,  the  residual  field  is  still  such  that  it 
makes  a  significant  difference  in  the  direct  integration  method  whether 
the  boundary  values  are  assumed  to  refer  to  the  geoid  or  to  some  mean 
elevation  in  the  area.  This  significant  sensitivity  to  data  height  error 
is  noted  for  altitude  below  30  km.  This  significant  sensitivity  is 
expected  since  the  field  roughness  has  not  been  removed  through  the 
subtraction  of  the  180*field.  The  field  roughness  directly  causes  high 
sensitivity  of  the  modeling  procedures  to  height  errors,  since  high 
frequency  field  irregularities  are  drastically  affected  by  upward  or 
downward  continuations.  Even  the  use  of  a  mean  elevation  surface  as 
assumed  reference  level  for  the  data  is,  of  course,  in  error  since  the 
anomaly  data  actually  refer  to  the  physical  topography  which  is  not  an 
equipotential  (level)  surface.  To  reduce  the  error  it  is  desirable  to  use 
the  residual  topographic  model  as  complementary  model  to  smooth  the 
field. 

Test  5  is  cf  interest  to  show  the  effects  of  residual  topographic 
masses  outside  the  2*x2*  area  surrounding  the  computation  points. 
These  outer  zone  effects  were,  on  the  average,  12%  of  the  total  residual 
topographic  effects  at  our  computation  points. 

Test  6  was  used  to  see  if  the  approximate  residual  topographic 
anomaly  ig^  of  (4.29)  can  be  used  in  place  of  the  rigorous  4g^  outside 
the  2*x2*  area  surrounding  the  computation  points.  The  resulting 
errors  from  this  test  can  be  considered  tolerable,  since  these  errors  are 
actually  the  effect  of  two  types  of  errors:  one  is  the  error  of 

approximating  Ag^  by  Ag^,  and  the  other  is  the  error  of  the  direct 
integration  method  in  representing  the  residual  topographic  field  (more 
about  this  in  Test  9).  It  is  important  to  be  able  to  approximate  Ag*=  by 
Ag‘  in  some  outer  zone  because,  for  example,  the  computation  of  rigorous 
values  of  Ag^  for  the  entire  7*x9*  test  area  would  lead  to  84x108=9072 
values  which  would  be  prohibitive  to  compute  using  Forsberg’s  prism 
integration  program.  The  computation  of  a  24  x  24  array  of  rigorous 
Ag‘  -values  inside  the  2*x2*  inner  zone  area  already  took  about  10 
minutes  of  CPU  time  on  the  OSU  Amdahl  system. 


Comparing  Test  7  with  Test  4  we  see  very  similar  errors,  supporting 
the  idea  that  most  of  the  errors  caused  by  varying  data  level  occur  in 
the  (high  frequency)  residual  topographic  part  of  the  total  field. 

Comparing  Test  8  with  Test  4,  we  see  significantly  smaller  errors  for 
Test  8.  This  is  due  to  the  smoothness  of  the  residual  integral  field 
used  in  Test  8.  The  low  sensitivity  of  the  residual  integral  field  with 
changing  data  level  means  that  now  the  assumption  of  the  integral 
model,  namely,  that  the  surface  data  refer  to  a  single  level,  will  produce 
only  small  errors.  It  is  indicated  from  Test  8  that  at  very  low  altitudes, 
below  10  km,  we  still  encounter  a  significant  although  small  model 
sensitivity  to  changing  data  level. 

Teat  9  shows  the  magnitudes  of  errors  that  can  be  encountered 
when  using  the  direct  integration  method  to  model  a  rough  field.  Part 
of  the  errors  can  be  attributed  to  the  inability  of  the  direct  integration 
method  to  properly  account  for  the  changing  level  of  the  given  anomaly 
data,  (i.e.,  to  properly  account  for  the  shape  of  the  topography). 

Note  that  in  Test  9  a  value  of  h„=1.5  km  was  used  in  (4.25).  In  a 
separate  testing  we  also  used  h,a=0,  and  generated  the  corresponding 
errors  as  in  Test  9.  This  was  done  to  confirm  that  the  use  of  the  mean 
elevation  hn^l.S  km  would  yield  smaller  errors  than  the  use  of  h^sO. 
The  statistics  of  the  errors  for  h4,=0  and  h«=1.5  km  are  shown  in  Table 
6.  We  see  a  clear  tendency  for  the  mean  and  standard  deviation  of  the 
errors  associated  with  hM=i.5  km  to  be  smaller  than  those  for  h,„-0  km, 
so  that  hai^l.S  km  is  the  preferable  value  to  use. 

We  close  this  section  by  giving  Table  7,  which  shows  the  final  model 
values  for  the  test  trajectory  through  point  A  (Figure  25).  The  total 
disturbance  values  were  obtained  by  adding  together  the  values  from 
the  spherical  harmonic,  residual  topographic,  and  residual  integral 
models.  Also  shown  in  Table  7  (row  V)  are  the  errors  of  the  direct 
integration  method  in  representing  the  residual  topographic  field,  as 
explained  in  the  description  for  Test  9.  To  get  the  disturbance  values 
for  the  case  when  the  residual  topographic  model  is  not  used  as  a 
complementary  model,  we  have  to  add  corresponding  numbers  given  in 
rows  IV  and  V  of  Table  7. 


Table  7.  Values  of  the  Disturbance  Components  in  Space  as  Evaluated 

from  the  Models  of  Chapter  4,  for  the  Test  Point  A  (4  =  33 'N, 
X  =  254’.5E).  Row  V  gives  the  errors  of  the  direct  integration 
method  in  representing  the  residual  topographic  model  (see 
Test  9).  Units:  mgals. 


MODEL 


(spherical 

hamonic 

model) 


(residual  .t 
topographic  * 

model)  .t 

«X 


(residual 
integral 
model ) 


(total) 


tf^)°  - 


D 

«x)  - 


5  km 

10  km 

30  km 

100  km 

26.31 

23.41 

14.56 

1.60 

8.52 

8.00 

6.41 

3.88 

-3.43 

-2.80 

-1.00 

1.14 

47.18 

39.62 

17.48 

-0.64 

-5.00 

-3.92 

-1.77 

-0.64 

-22.50 

-14.59 

-1.73 

0.48 

-12.50 

-7.72 

-0.71 

2.60 

4.97 

3.46 

2.01 

1.31 

-4.97 

-2.52 

-2.09 

-1.64 

60.99 

55.31 

31.33 

3.56 

8.49 

7.54 

6.65 

4.55 

-30.90 

-19.91 

-4.82 

-0.02 

0.35 

-0.68 

-1.25 

-0.62 

-1.28 

0.01 

0.19 

-0.05 

8.18 

4.85 

1.83 

0.34 

1.48 


5.  DISCRETB  APPROACHES  TO  MODELING  THE  SPATIAL  DISTURBANCE 

VECTOR 

Discrete  approaches  to  modeling  the  external  field  are  based 

on  solving  the  so-called  Bjerhammar  problem.  Bjerhammar  (1964) 
formulated  the  modeling  problem  in  this  way:  given  a  finite  number  of 
gravity  anomalies  at  points  of  a  non-spherical  surfacsi  find  an  analytical 
approximation  to  the  disturbing  potential  T  in  space  such  that  the  given 
anomaly  vidues  are  satisfied.  Solutions  that  have  been  developed  for 
the  Bjerhammar  problem  now  include  possibilities  to: 


1.  satisfy  heterogeneous  data  types  (e.g.,  gravity  anomalies,  deflections 
of  the  vertical,  gravity  gradients)  at  the  given  observation  points; 

2.  filter  the  data  to  minimize  the  influence  of  measurement  errors; 

3.  remove  systematic  influences  that  make  T  non-harmonic,  such  as 
time  dependent  variations,  reference  system  errors  in  orientation 
and  geocentricity,  and  station  position  errors  (Tscherning,  1981). 


Well-known  solutions  presented  for  the  Bjerhammar  problem  are 
based  on  collocation,  which  is  a  numerical  method  for  fitting  an 
analytical  approximation  of  a  function  to  a  finite  number  of  given 
functionals  of  the  function  (see  Moritz,  1978b;  see  also  Moritz,  1978a,  p. 
32,  for  a  simple  example  of  collocation).  In  fact,  the  Bjerhammar  problem 
as  stated  above  reads  as  a  collocation  problem.  Details  of  the 
application  of  collocation  in  solving  the  Bjerhammar  problem  have  been 
presented  under  different  analytical  structures.  Original  presentations 
are  given  by  Krarup  (1969)  using  Hilbert  spaces,  Moritz  (1970a)  using 
Wiener  prediction  and  filtering,  and  Bjerhammar  (1975)  using  Dirac 
functions.  Extensive  treatment  of  various  aspects,  both  for  practical 
application  and  fcr  deeper  analytical  understanding  of  collocation,  are 
given  by  Moritz  (1980).  Review  papers  dealing  with  the  theory  of 
collocation  are  given  by  Moritz  (1978a,  1978b)  and  Bjerhammar  (1978). 
Tscherning  (1981)  gives  a  review  paper  comparing  the  theory  and 
computational  properties  of  different  gravity  field  representations, 
including  collocation  and  its  combination  with  other  methods,  and 
summarizing  also  the  results  of  previous  numerical  investigations  by 
various  authors. 

Of  special  interest  are  the  numerical  investigations  of  Sjoberg 
(1975,  1978),  who  used  both  model  computations  and  limited  real  data 
computations  to  predict  gravity  anomalies  and  deflections  of  the  vertical 


on  the  surface  of  the  topography,  starting  from  given  surface  gravity 
anomalies.  He  concludes  that  similar  prediction  results  are  obtained 
from  the  use  of  different  base  functions  in  the  collocation  procedure 
(details  of  the  procedure  are  given  below),  provided  that  the  radius  of 
the  internal  sphere  is  properly  chosen  for  each  base  function.  In  this 
chapter  we  will  investigate  the  validity  of  this  conclusion  for  gravity 
disturbance  vector  computations  in  space. 

Of  interest  to  us  also  are  the  numerical  investigations  of 
Katsambalos  (1981)  who  performed  simulation  studies,  precisely  on  the 
computation  of  the  gravity  disturbance  vector  in  space  from  surface 
gravity  anomaly  data.  He  used  model  values  using  the  "Molodensky 
mountain  model"  (Molodensky  et  al.,  1962).  Key  conclusions  of 
Katsambalos  were  questioned  by  Tscherning  (1983a)  with  regards  to  the 
least  squares  collocation  approach,  and  by  Bjerhammar  and  Sjoberg 
(1982,  private  communications)  with  regards  to  the  so-called  Dirac 
approach.  Both  the  least  squares  collocation  and  Dirac  approaches  are 
collocation  techniques,  as  we  will  detail  below.  The  results  of  our  own 
numerical  experiments,  to  be  presented  in  this  chapter,  support  the 
validity  of  the  questions  raised.  Specifically,  our  results  show  the 

following: 

1.  the  least  squares  collocation  and  Dirac  approaches  can  give  about 

the  same  quality  of  results  if  the  free  parameters  of  the  base 

function  used  (degree  variances,  radius  of  the  Bjerhammar  sphere) 
are  properly  selected. 

2.  the  conclusions  of  Katsambalos  based  on  his  Tables  9.3  and  9.4 

(ibid.,  p.  116),  namely,  that  the  "Dirac  approach  should  not  be  used 
for  the  interpolation  of  gravity  anomalies  on  the  surface  of  the 
earth,  unless  the  data  is  sufficiently  dense  to  ensure  an  accurate 
prediction",  should  be  modified.  His  Table  9.3  was  incorrectly 
computed,  with  the  correct  table  being  given  in  Sjoberg  (1978,  p. 
68,  Table  8.3).  Kis  Table  9.4  represents  only  a  very  specific 

example,  and  cannot  be  generalized  to  form  the  conclusion  that  the 
Dirac  approach  tends  to  predict  zero  values  in  between  the  data 
points.  The  conclusion  should  rather  be,  that  the  Dirac  approach 
can  be  used  for  interpolation  provided  the  radius  of  the  internal 
sphere  is  properly  selected.  The  large  prediction  errors  noted  by 
Katsambalos  (ibid.,  p.  113)  when  doing  real  data  computations  over 
Canada  was  due  to  the  relatively  large  spacing  between  the  given 
data  points  (see  Figure  28  of  this  chapter),  and  not  due  to  a 
deficiency  in  the  least  squares  collocation  or  Dirac  prediction 
methods  as  Katsambalos  suggested.  We  will  present  additional 
computations  using  the  mentioned  real  data  in  Canada. 


5.1  The  Collocation  Procedure 

We  already  mentioned  that  the  word  collocation  denotes  the 
construction  of  an  approximation  to  a  function  by  fitting  an  analytical 
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approximation  of  the  function  to  a  finite  number  of  functionals  of 

the  function  (see  Moritz,  1978a,  1978b,  for  example).  In  gravimetric 
geodesy  the  fundamental  function  to  be  approximated  is  usually  the 
disturbing  potential  T,  and  the  given  functionals  are  (for  example) 
gravity  anomalies  Ag,  deflections  of  the  vertical  (  and  rj,  and  height 
anomalies  f,  which  are  all  related  to  T  through  well-known  linearized 
relations  in  spherical  approximation: 


(5.1) 


f  =  ■ 


JL  ll 

ry 


(5.2) 


1  »T 
ry  cos  ▼ 


(5.3) 


To  briefly  explain  the  use  of  collocation  in  constructing  an 
approximation  to  the  disturbing  potential  T  in  space,  we  use  the 
presentation  of  Moritz  (1978b).  It  is  first  postulated  that  T  can  be 
approximated  by  a  function  f  in  the  form: 


T(P)  «  f(P)  =  E  bj  ♦j  (P). 
j=i 


(5.4) 


where 

P  computation  point  in  space 

bj  constant  collocation  coefficients  to  be  determined 
♦j  a  priori  selected  base  functions 

It  is  considered  that  T  is  harmonic  outside  the  earth's  masses,  and  to 
enforce  harmonicity  on  the  approximating  function  f,  the  selected  base 
functions  are  harmonic: 

=  0  ,  (5.5) 

where  A  denotes  here  the  Laplace  operator.  For  the  determination  of 
the  collocation  coefficients  bj  it  is  assumed  that  errorless  values  of 
are  given  at  k  spatial  points  on  or  above  the  earth’s  surface,  where  the 
1 1  are  linear  functionals  (Moritz,  1980,  p.  36)  of  T  or  f: 
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1)  =  L)T  =  L(f  , 


(5.6) 


where  L  (  denotes  the  linear  operation  to  transform  T  (or  f )  to  <  j .  For 
example,  (5.1)  is  in  the  form  (5.6)  with 


(5.7) 


=  -  >  -1 
l«r,  r,J 


(5.8) 


The  collocation  condition  is  now  the  exact  reproduction  of  the  given 
by  the  function  f,  i.e.,  substituting  (5.4)  into  (5.6): 


r~  Lt  ♦j  (Pt)  =  . 

j=i 


(5.9) 


In  matrix  notation,  condition  (5.9)  becomes: 


B  b  =  4  , 


where  the  elements  of  the  B-matrix  are: 


B, i  =  L<  ♦,  (Pi)  . 


(5.10) 


(5.11) 


i.e.,  the  operator  Li  is  applied  to  the  base  function  4;  and  the  result 
evaluated  at  the  point  P|  where  *i  is  given.  From  (o.lO)  we  formally 
get  the  collocation  coefficients: 


b  =  B"‘l  . 


(5.12) 


In  actual  implementation,  (5.10)  may  be  solved  for  b  by  iteration  or  by  a 
linear  equation  solution,  instead  of  the  direct  matrix  inversion  in  (5.12). 

Knowing  the  collocation  coefficients  bj,  (5.4)  gives  the  required 
approximation  to  T.  Any  other  functional  LT  of  T  can  be  found  in  space 
by  applying  the  operator  L  to  (5.4): 


95 


LT(P)  s  Lf(P)  =  ^  bj  L*j  (P) 
j=l 


(5.13) 


For  example,  in  computing  the  components  of  the  gravity  disturbance 
vector  in  space  we  have  the  linear  relations: 


-  -il 

-  77 


.  1  *T 


1  ^ 
r  cos  ? 


(5.14) 


(5.15) 


(5.16) 


Note  that  (5.13)  is  also  of  the  form  (5.4),  with  the  fundamental  function 
T  replaced  by  the  derived  function  LT,  and  the  fundamental  base 
function  ♦j(P)  replaced  by  the  "propagated"  base  function  L*j(P).  This 
completes  our  brief  account  of  the  use  of  collocation  in  approximating  T 
and  quantitites  derived  from  T  (see  (5.4)  and  (5.13)). 

An  important  issue  in  the  above  collocation  procedure  is  the  choice 
of  base  ,  function.  In  this  report  we  will  examine  two  major  approaches 
in  selecting  the  base  function,  namely,  the  Dirac  approach  (Bjerhammar, 
1975;  Sjoberg,  1978)  and  the  least  squares  collocation  approach  (see, 
e.g.,  Moritz,  1978b).  We  realize  from  the  geodetic  literature  on 
collocation  that  a  set  of  equations  purportedly  defining  one  method  can 
in  certain  cases  be  interpreted  as  defining  another  method.  Examples  of 
this  are  found  in  Sjoberg  (1978,  pp.  9,  12),  Tscherning  (1983a), 

Lelgemann  (1981),  Brennecke  and  Lelgemann  (1984),  Bjerhammar  (1975, 
1978).  For  our  purposes  we  will  distinguish  between  the  Dirac  approach 
and  least  squares  collocation  approach  in  a  simple  way,  by  specifying 
the  equations  we  are  using  under  each  method.  Then  we  will  conduct 
numerical  experiments  using  the  equations  and  compare  the  results.  The 
reader  interested  in  a  theoretically  refined  distinction  between  the  Dirac 
and  least  squares  collocation  approaches,  and  in  their  relationship  to 
each  other  and  to  other  models,  should  consult  the  references  cited  and 
Moritz  (1980,  starting  from  Sections  12,  30,  38). 


5.2  The  Dirac  Approach  to  Collocation 

The  Dirac  approach  to  collocation  has  been  interpreted  as  a 
generalized  point  mass  modeling  (Sjoberg,  1978,  p.  13;  see  also  Forsberg, 
1984,  pp.  24-27;  and  Brennecke  and  Lelgemann,  1984).  In  this 
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with 


cos  V’pj  =  sin  ♦p  sin  ♦j  +  cos  ♦p  cos  ♦j  cos  (Xj  -  Xp)  ,  (5.19) 

4,  X  ...  geodetic  latitude  and  longitude. 

The  base  function  (5.18)  is  in  the  form  of  the  kernel  function  E  used  in 
the  Brovar  solution  mentioned  earlier  (see  Moritz,  1980,  p.  366). 

Various  choices  of  the  k„  in  (5.18)  result  in  various  further 
interpretations  of  the  quantities  Uj*  appearing  in  (5.17).  To  see  this  we 
first  create  from  the  impulses  Uj*  a  function  u*  on  the  Bjerhammar 
sphere  (Bjerhammar,  1975,  (3)): 


k 

u*(Q)  =  4nRB»  ^  Uj*  «(Q-j)  ,  (5.20) 

j=l 

where 

0  arbitrary  point  on  the  Bjerhannar  sphere 

5(0-j)  Dirac  delta  function  (see,  e.g.,  Moritz,  1980,  pp. 

27-30).  This  function  has  the  characteristics  that 
tf(Q-j)=0  if  Q^j,  and  that  the  integral  of  tf(O-j) 
is  equal  to  one  for  every  neighborhood  of  point 
j  howf  ver  saall  this  neighborhood  is  chosen  to  be. 


Then,  (5.17)  can  be  rigorously  written  in  the  integral  form  (compare 
Bjerhammar,  1975,  (l)-(5a)  or  Katsambalos,  1981,  (5.33)-(5.35): 


T(P)  *  f(P)  =  ^  f  J  u*(Q)  base(T(P),  u*(Q))  da(Q)  ,  (5.21) 

a 


where  the  integration  is  carried  over  the  unit  sphere  <r  and  the  integral 
kernel  is,  from  (5.18): 


base(T(P),  u*(Q))  =  Rg  K  Pn(cos  ipg)  .  (5.22) 

n=2  ’’ 


From  the  theory  of  eigenvalues  and  eigenfunctions  of  isotropic  integral 
operators  (see,  e.g.,  Jekeli,  1980)  we  can  write  from  (5.21)  and  (5.22)  the 


following  relation  between  the  n^"  surface  harmonics  of  T*  and  u*  on 
the  Bjerhammar  sphere: 


T„*  =  k„  u„«  .  (5.23) 

(The  denotes  reference  to  the  Bjerhammar  sphere).  The  last  equa¬ 
tion  reveals  directly  what  type  of  qtiantity  the  u*  is,  in  relation  to  the 
disturbing  potential  T.  For  example,  if  k„=l  is  chosen  then  the  spectral 
relation  (5.23)  reveals  that  the  u*  being  solved  on  the  Bjerhammar 
sphere  is  the  gravity  anomaly  (see  (5.24)  below);  from  (5.20)  the  Uj*  in 
this  case  can  be  termed  "gravity  anomaly  impulses."  We  have  from 
(5.23): 


=  k„  u„*  .  (5.23a) 

In  this  report  we  will  numerically  experiment  on  three  types  of 
choices  for  the  k„  in  (5.22).  These  choices  will  be  those  that  imply  the 
solution  on  the  Bjerhammar  sphere  of  the  gravity  anomaly  ^g*,  the 
single-layer  density  and  the  double-layer  density  Ji*.  More  details 
on  these  quantities  follow: 

1.  the  gravity  anomaly  ta*.  Here,  the  Agj^«  are  gravity  anomaly 
impulses.  This  is  the  original  system  used  by  Bjerhammar  (1964),  and 
for  which  numerical  experiments  were  performed  by  Sjoberg  (1975,  pp. 
86-89;  1978,  pp.  64-69)  and  Katsambalos  (1981)  under  the  label  of  "Dirac 
approach."  From  Heiskemen  and  Moritz  (1967,  p.  89)  we  have  the 
spectral  relation  between  T*  and  dg*: 


T„* 


(5.24) 


Conparing  (5.23)  and  (5.24)  the  choice  of  k„  in  this  case  is: 


k„(dg*)  =  1  . 


(5.25) 


The  base  function  (5.18)  then  becomes  the  extended  Stokes  function, 
with  closed  form  (see  ibid.,  (6-30),  (6-45): 


base(T(P),  Agj*)  =  Rgt  [|  1-3D-  t  cos^  (5  +  3ln  ^  t  cos^^pj  j 


where 


(5.26) 
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Tp  =  R+hp  (R  =  6371  km;  hp  =  height  of  P  above  the  geoid)  (5.28) 


cos  i)  =  cos  ^pj 


(see  (5.19)) 


1  -  2t  cos  i)  +  . 


(5.29) 


(5.30) 


2.  the  single-layer  density  u*.  Here,  the  pj*  are  point  masses. 
This  system  (of  Dirac  approach  to  collocation)  is  discussed  in  Sunkel 
(1981b,  1983b).  Prom  Sunkel  (1981b,  (2.5))  we  get  the  following  spectral 
relation  between  T*  and  p*: 


^  RB(2n+l)  • 


(5.31) 


Compeuring  (5.23)  and  (5.31)  the  choice  of  k„  in  this  cgfe  is: 


k  (u*)  =  — — 
Rg»(2n+1) 


(5.32) 


The  base  function  (5.18)  becomes  simply  the  reciprocal  distance  function, 
with  closed  expression  (Heiskanen  and  Moritz,  1967,  p.  35): 


baae(T(P),  pj*)  =  ^  ^  P„(cos  i)p^)  =  -j 


ia  i--  ^rp^ 
n=0 


(5.33) 


where 


rp  Rg  -  2  Rg  rp  cos  Vpj 


(5.34) 


Note  that  for  exact  agreement  of  (5.18)  and  (5.33)  we  should  start  the 
summation  in  (5.18)  from  n=0. 

3.  the  double- layer  density  u».  Here,  the  lij*  are  point  dipoles. 
The  concept  of  point  dipoles,  and  in  general  multipoles,  is  discussed  in 
Meissl  (1981,  pp.  184-190).  The  base  function  corresponding  to  the  use 
of  multipoles  on  the  Bjerhammar  sphere  can  be  expressed  as  a 
derivative  of  (5.33)  with  respect  to  Rg  (ibid.,  7.7)): 


-r  [(i)]  , 

n=0 


(5.35) 


w-g-w  y-w  imr*^  j 


■jr^ jr».'c^.\’^  y  ■'_’»' 


J■."'.’v^.■v^^Jr-  .r.  irwTrvirvTr^ 
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where  the  order  N  of  differentiation  is  the  order  of  the  multipole,  e.g., 
N=1  for  a  dipole  and  N=0  for  a  monopole  (i.e.,  a  point  mass).  Using 
(5.35)  with  N=l,  we  get  the  following  base  function  for  our  dipole 
system,  in  series  form: 


base(T(P),  >Ij*)  =  )  n  [^]  P„(co8  Vpj)  , 

n=l 


(5.36) 


and  in  closed  form  (differentiating  (5.34)); 


base(T(P),  Mj#)  =  -«* 


(5.37) 


Comparing  (5.18)  and  (5.36)  gives  the  implied  k„: 


k  (u*)  =  — 

^  RB’(2n+l) 


(5.38) 


Using  (5.38)  into  (5.23)  we  get  the  spectral  relation  between  the 
disturbing  potential  T*  and  the  double-density  layer  Jl*  on  the 
Bjerhammar  sphere: 


T„0  = 


RB*(2n+l) 


Mn*  • 


(5.39) 


We  have  now  completed  the  description  of  the  various  Dirac  systems 
that  we  will  use  in  our  numerical  experiments.  These  systems  have  been 
defined  through  particular  choices  of  k„  in  the  base  function  (5.18).  In 
principle,  the  choice  of  k„  in  (5.18)  is  arbitrary,  provided  only  that 
(5.18)  converges  and  that  we  can  solve  for  the  collocation  cefficients  bj 
in  the  collocation  formula  (5.4)  (details  of  the  computation  of  the  bj  are 
given  later).  We  have  chosen  the  k„ -expressions  (5.25),  (5.32),  and 
(5.38),  leading  respectively  to  the  gravity  anomaly  impulses,  point 
masses  and  point  dipole  systems,  because  the  said  systems  are  relatively 
well-known  in  the  geodetic  literature.  More  numerical  experiments  are 
needed  on  these  systems,  especially  using  real  data  and  computations  in 
space. 


Note  that  both  the  k„(4g«)  of  (5.25)  and  k„(M*)  of  (5.32)  are 
proportional  to  a  constant  for  n  whereas  the  kn(?*)  of  (5.38)  is 

proportional  to  n  for  n  -»  *.  Adopting  the  terminology  of  Lelgemann 
(1981),  the  order  of  k„(&g*)  and  are  equal  to  zero,  whereas  the 

order  of  kn(JZ«)  is  equ^  to  1.  The  order  of  k„  in  n,  denoted  by  0(k„), 


aS 


I 


■i  *  to  ' 

;i5l 


I 


is  defined  as  "the  highest  power  of  n  within  the  explicit  expression  for 
kn»  e.g.,  for  kn  =  (2n-fl)/(n-f  1)  we  obtain  0(k„)  =  0"  (ibid.).  Analogous  to 
the  empirical  findings  of  Lelgemann  in  his  interpolation  problems,  we 
also  found  that  the  0(kn)  strongly  influences  the  optimal  radius  Rg  that 
should  be  chosen  for  good  collocation  results,  and  that  systems  with  the 
same  0(k„)  have  about  the  same  optimal  Rg.  The  Ag*-  and  ;i*-systems 
required  about  the  same  depth  to  the  Bjerhammar  sphere,  these  two 
systems  having  the  same  0(k„).  The  Ji*-system  required  a  deeper 
Bjerhammar  sphere,  because  this  system  had  a  larger  0(k„)  than  the 
other  two  systems.  From  (5.23a)  one  can  see  that  the  larger  the  0(k„) 
the  smoother  (less  high  frequency  energies)  the  u<  relative  to  Ag*,  with 
0(k„)=0  implying  the  same  smoothness  for  u*  and  Ag«.  One  should 
remember  that  spectral  multiplication  by  n  implies  amplification  of  high 
frequency  energies  of  a  function. 

Let  us  now  turn  to  the  computation  of  functionals  of  T,  from  the 
general  representation  (5.17)  of  the  Dirac  approach.  In  accordance  with 
(5.13),  a  functional  LT  of  T  can  be  found  by  applying  the  operator  L  to 
(5.17): 


k 

LT(P)  *  Lf(P)  =  ^  base(T(P),  Uj*)  .  (5.40) 

j=l 


FroB  (5.40)  we  can  write: 


base(LT(P),  Uj*)  =  L  base(T(P),  Uj») 


(5.41) 


The  last  equation  expresses  a  "propagation  of  base  function,"  analogous 
to  the  propagation  of  covariance  functions  well-known  in  least  squares 
collocation  (see  Moritz,  1980,  p.  86).  As  an  example  of  (5.40)  and  (5.41) 
let  LT(P)  =  Ag(P)  and  Uj*  =  Agj*;  using  (5.40),  (5.41),  and  the  operator 
(5.8)  for  Ag  we  have: 


k 

4g(P)  =  base(Ag(P),  Agj*)  , 

j=l 


(5.42) 


where  the  propagated  base  function  is: 


base(Ag(P),  Agj*)  =  ■—  -  base(T(P),  Agj*) 


(5.43) 
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and  the  fundamental  base  function  base(T(P),  ^gj*)  is  given  by  (5.26). 
In  the  following  we  show  the  propagated  base  functions  (5.41)  for  our 
gravimetric  quantities  of  interest  (i.e.,  for  LT=Ag,  fip,  fiy,  and  6\),  and 
under  our  various  Uj*-9ystems  (i.e.,  UjtzAg^t,  /ij*,  Jij»). 

1.  igj*  -system.  Formally  we  can  apply  (5.41)  to  get  the  propagated 
base  functions  for  this  system.  However,  we  already  know  these 
functions  to  be  the  Poisson  kernel  for  the  upward  continuation  of 
gravity  anomalies  (Heiskanen  and  Moritz,  1967,  sec.  2-15),  and  the 
kernels  for  the  "direct  method"  of  horizontal  and  radial  disturbance 
computations  (ibid.,  sec.  6-4).  We  have: 

baae(4g(P),  Agj*)  =  ^  ^  ;  (5.44) 

baae(Ap(P),  Agj»)  =  |  +  1  -  6D  - 

-t  cos  ^  (l3  ^  6  ID  1  -  t..cos  .  Dj]  .  (5  45) 


base(4^(P),  Agj*) 

cos  Ot 

base(Ax(P) »  Hj*) 

sin  a 

t’  sin 


vlsl 


.1-8-3  1  -  t  °  -  3  1  t  ° 

D  D  sin*  V*  2 


(5.46) 


In  (5.46)  a  is  the  clockwise  azimuth  from  North,  from  the  projection  of  P 
on  the  Bjerhammar  sphere  to  the  point  of  Agj*.  We  have  (Rapp,  1966b, 

p.  17): 


cos  a  =  cos 


cos  ♦p  sin  ♦j  ~  sin  cos  ♦j  cos(Xj  -  X, 

sin  V' 


cos  sin(.\i  -  X, 

sin  <*  =  sin  ap,  =  - '  i  ‘ 

sin  if 


(5.47) 

(5.48) 


2.  jijS-system  (Point  Masses).  For  this  system  we  need  to  apply 
the  operators  in  (5.1),  (5.14),  (5.15),  (5.16)  to  the  fundamental  base 


function  (5.33)  to  obtain  the  propagated  base  functions.  The  results 


ba,e(Ae(p).  PJ.)  =  ('  ~  -  i)  ; 

bMe(«,(P).  „J»)  .  ^  [1  -  t^CMyj  , 


base(tfg(P),  fij*) 


baae(6x(P).  Mj*) 


t  simfr 
rp»  D’ 


(4.49) 


(5.50) 


(5.51) 


In  deriving  (5.51)  we  have  used  (Heiskanen  and  Moritz,  1967,  p.  234): 


=  -c 


^  =  -c 

*X 


(5.52) 


3.  /ijt-system  (Point  Dipoles).  Again  we  apply  the  operators  in 
(5.1),  (5.14),  (5.15),  (5.16),  this  time  to  the  fundamental  base  function 
(5.37)  to  obtain  the  propagated  base  functions  for  this  system.  The 
results  can  be  written  as  follows: 


b..e(4j(P),  iij.)  .  -L  cos^f  -  CO.  f  _  2.^ 


1  (3  cos  9  cos  B  -  cos 


p-^)  ;  (5.53) 


base(fi,.(P),  Mj*)  =  ^  ( 


base(6f(P),  JIj*) 


base(fix(P).  Mj*) 


(5.54) 


^  (3_S_p^  ^  ^ 


Since  the  fundamental  base  function  (5.37)  for  the  point  dipole  system 
has  been  found  by  differentiating  the  base  function  (5.33)  for  the  point 
mass  system  with  respect  to  the  radius  Rbi  then  it  follows  from  the  idea 
of  base  function  propagation  that  (5.53)  -  (5.55)  can  also  be  obtained  by 
differentiating  (5.49)  -  (5.51)  with  respect  to  Rg. 

Finally,  to  complete  our  description  of  the  Dirac  approach  to  the 
collocation  problem,  we  need  to  discuss  the  computation  of  the 
coefficients  Uj<  (the  generalized  point  masses)  appearing  in  (5.17).  We 
will  determine  the  u.*  (j  =  1,  2,  ...,  k)  using  k  given  functionals  LT  of 
T  on  the  earth’s  surface,  and  an  inversion  of  the  k  x  k  system  of  linear 
equations  arising  from  (5.40).  In  the  principle  the  Dirac  system  can 
handle  heterogeneous  data  types,  i.e.,  different  types  of  functionals  of  T 
such  as  gravity  anomalies,  deflections  of  the  vertical,  and  height 
anomedies.  In  principle  also,  a  filtering  of  observational  errors  can  be 
built  into  the  method  by  using  more  than  k  given  functionals  and 
performing  a  least  squares  fit  to  the  data  based  on  (5.40). 

However,  our  study  will  be  limited  to  the  use  of  a  homogeneous  data 
type,  the  gravity  anomaly,  and  a  unique  solution  of  k  x  k  system  of 
linear  equations  to  solve  for  the  Uj*.  The  gravity  anomaly  is  still  the 
most  common  gravimetric  data  type,  being  the  observable  arising  from 
the  Molodensky  formulation  (Molodenakii,  ei  al.,  1962)  of  the  geodetic 
boundary  value  problem  (g.b.v.p.).  Bjerhammar  and  Svenson  (1983) 
discuss  the  straightforward  application  of  the  Dirac  approach  to  the 
so-called  "fixed  boundary"  formulation  of  the  g.b.v.p.,  where  the 
observables  are  gravity  disturbances.  Regarding  the  use  of  a  unique 
solution  for  the  uj*  instead  of  a  least  squares  solution,  we  believe  that 
it  would  be  best  to  separate  the  problem  of  modeling  the  external  field 
from  surface  data  on  the  one  hand,  and  the  problem  of  the 
establishment  of  an  optimal  set  of  surface  data  on  the  other  hand.  The 
latter  problem  can  be  handled  in  an  a  priori  step  involving  editing, 
prediction,  filtering,  and  adjustment  of  observational  data. 

To  solve  for  the  collocation  coefficients  Uj*  (j  =  1,  2,  ...,  k)  from  k 
given  gravity  anomaly  values  on  the  surface  of  the  earth,  we  first  form 
the  k  X  k  system  of  lineeu*  equations  (see  (5.40),  (5.41)): 

k 

Agj  =  y  Uj*  base(Agi,  Uj»),  i  =  1,  2,  k  ,  (5.58) 

j=l 


where  in  our  case,  base  (Ag,,  Uj*)  can  be  any  of  (5.44),  (5.49),  or  (5.53). 
Since  the  system  (5.58)  tends  to  be  large,  we  employ  an  iterative 
solution.  The  number  of  iterations  required  for  a  convergent  solution 
indicates  the  degree  of  stability  of  the  equation  system.  We  will  use  the 
Gauss-Seidel  iteration  method  employed  also  by  Katsambalos  (1981,  p.  64). 
Letting  (m)  denote  the  current  iteration  step,  we  have: 


A, j  =  base(Ag, ,  Uj*)  . 


(5.60) 


Note  from  (5.59)  that  the  most  recently  computed  value  of  Uj«  is  used  at 
each  iteration  step.  The  gravity  anomalies  implied  by  a  current 
iteration  step  are: 

k 

Ag/")  A,j  (i  =  1,  2,  ....  k)  ,  (5.61) 

j=l 

and  the  current  residuals  are: 

*,(<")  =  Ag,  -  Ag,("')  ,  (i  =  1,  2,  ....  k)  .  (5.62) 

It  is  reasonable  to  terminate  the  iterations  when  at  least  one  of  these 
conditions  is  met  (Sjoberg,  1978t  p.  68): 
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®Rxs  ; 

Max 
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(5.63) 
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where  the  tolerances  crhs  (RMS  residual),  (maximum  residual),  and 

m,aaK  (inaximum  number  of  iterations)  are  to  be  prescribed,  for  example 
(ibid.): 


Cfiifs  —  0.25  Bgal 
=  0.5  ngal 


(5.64) 


This  completes  our  description  of  the  Dirac  approach  to  collocation. 
To  summarize,  the  formal  procedures  in  the  Dirac  approach  is  to  model 
T(P)  as  the  sum  of  potentials  of  k  generalized  point  masses  (see  (5.17)) 
that  are  distributed  on  an  internal  sphere  of  chosen  depth  (a 
Bjerhammar  sphere  of  radius  Ro).  The  values  of  the  "point  masses"  can 
be  solved  from  k  given  functionals  LjT(Pj)  of  T  on  the  earth's  surface 
by  an  inversion  of  a  k  x  k  linear  equation  system  arising  from  (5.40)  in 
general,  or  (5.58)  in  our  particular  case.  Various  choices  of  the  free 
factor  k„  in  the  base  function  (5.18)  lead  to  different  types  of  "point 
masses"  on  the  Bjerheunmar  sphere.  We  decided  to  use  gravity  anomaly 
impulses,  the  usual  point  masses,  and  point  dipoles  as  generalized  point 
masses  for  our  later  numerical  experiments. 

In  the  above  Dirac  approach,  and  in  the  least  squares  collocation 
approach  to  be  discussed  in  the  next  section,  it  is  important  that  we 
also  use  the  spherical  harmonic  and  residual  topographic  models  of 
Chapter  4  as  complementary  models.  The  reason  is  the  limitation  in  the 
size  k  of  the  linear  system  (5.58)  that  can  be  feasibly  solved.  Since  k 
is  equal  to  the  number  of  data  points,  then  it  is  important  to  (a)  use 
the  spherical  harmonic  model  to  account  for  remote  zone  effects,  thereby 
reducing  the  required  data  coverage  for  the  collocation  solution;  and  (b) 
use  the  residual  topographic  model  to  account  for  very  detailed  gravity 
field  information,  thereby  reducing  the  required  density  of  data  for  the 
collocation  solution. 


5.3  The  Least  Squares  Collocation  Approach 

In  the  Dirac  approach  to  the  collocation  problem,  the  base  functions 
used  are  the  k  harmonic  functions  ba8e(T(P),  Uj*)  relating  the 
disturbing  potential  T(P)  to  each  of  the  k  Dirac  impulses  Uj*  on  the 
Bjerhammar  sphere  (see  (5.17)).  In  the  least  squares  collocation 
approach,  the  base  functions  to  use  are  the  k  empirical  covariance 
functions  between  T(P)  and  each  of  the  k  given  functionals  LjT(J)  on 
the  earth’s  surface.  This  choice  of  base  functions  arises  from  adding  to 
the  collocation  condition  (5.9)  the  condition  that  the  RMS  difference 
between  the  approximation  function  f  and  the  true  function  T  be 
statistically  minimized  (see,  e.g.,  Moritz,  1978a,  1978b).  Apart  from  this 
statistical  interpretation,  the  least  squares  collocation  approach  can  also 
be  justified  as  a  purely  analytical  approximation  technique,  in  the 
context  of  the  more  general  "minimum  norm  collocation"  in  a  Hilbert 
space  of  harmonic  functions  (see,  e.g.,  Tscherning,  1977)  or  in  the 
context  of  the  minimum  norm  Bjerhammar  solution  discussed  by  Sjoberg 
(1975,  1978).  Our  interest  in  this  study  is  to  first  give  the  equations 
defining  our  use  of  the  least  squares  collocation  approach,  then  later 
perform  numerical  experiments  using  these  equations. 

We  have  the  following  representation  in  the  least  squares 
collocation  approach: 


T(P)  =  f(P)  =  r  bj  cov((T(P).  LjT(J))  , 

j=l 


(5.65) 


where 


bj  collocation  coefficients 


cov(T(P),  LjT(J)) 


empirical  covariance  between  T(P)  and  the 
given  functional  LjT  at  point  J.  This  is 
derived  from  an  analytical  expression  of 
suitable  form  (see  (5.71))  fitted  to  the 
covariance  behavior  of  real  data  in  the 
area  of  computation. 


As  already  mentioned,  the  given  functionals  in  our  case  are  homo¬ 
geneously  the  gravity  anomalies  Agj  on  the  earth’s  surface  so  that  we 
wiU  use  (5.65)  in  the  form: 


T(P)  »  f(P)  =  E  bj  cov(T(P),  Agj) 
j=l 


(5.66) 


Any  functional  LT  of  T  can  be  expressed  by  applying  the  operator  L  to 
(5.66)  and  employing  covariance  propagation  (Moritz,  1980,  p.  86): 


LT(P)  s  Lf(P)  =  E  bj  cov(LT(P),  Agj)  . 
j=l 


(5.66a) 


For  the  sake  of  our  discussions  the  covariance  function  in  (5.66)  is  now 
written  in  the  general  form: 


cov(T(P),  Agj)  =  ^  d^  Ir^l"  ^  Pn(cos  V-pj)  ,  (5.67) 


where 


ro  chosen  radius  of  the  internal  sphere  (BJerhammar  sphere) 


v;v>vv;'>,fibv 


non-negative  factors  to  be  chosen,  function  of  n- 


d„ 


The  covariance  function  (5.67)  is  in  the  same  form  as  the  base  function 
(5.18)  of  the  Dirac  approach,  except  for  the  replacement: 

Rb  ^  •  (5.68) 


The  implied  auto-covariance  function  of  Ag  follows  from  (5.67)  using 
covariance  propagation  (ibid.)  and  the  spectral  relation  between  T  and 
Ag  (see  (5.24): 


cov(Ag(P),  Agj)  =  J  (2n+l)  P„(cos  V-pj)  .  (5.69) 

n=2  ■* 


Similarly,  we  obtain  the  implied  auto-covariance  function  of  T: 


cov(T(P),  Tj)  =  d*  Pn(cos  Vpj)  .  (5.70) 

n=2  ^  ^ 

The  last  equation  is  in  the  proper  general  form  of  the  homogeneous  and 
isotropic  auto-covariance  function  of  the  disturbing  potential  (see,  e.g., 
Tscherning  and  Rapp,  1974): 


cov(T(P),  Tj)  =  C- 
n=2 


f  rn^n-*-! 
Irprji 


P„(cos  V-Pj) 


(5.71) 


where  the  are  the  degree  variances  of  T  on  the  sphere  of  radius  ro. 
From  (5.70)  and  (5.71): 


<^n 


dn" 


2n-t-l 

(n-l)s 


(5.71a) 


For  local  applications  of  least  squares  collocation,  the  analytical 
expression  being  used  for  the  covariance  function  should  be  fitted  to 
the  local  covariance  behavior  of  the  data  in  the  area  of  computation. 
Usually,  a  fit  is  made  to  Ag-data,  these  being  the  most  readily  available 


data.  In  this  case  the  free  parameters  d„  and  ro  should  be  adjusted 
together  so  that  (5.69)  will  have  the  desired  fit  to  the  data.  Moritz 
(1980,  pp.  174-177)  has  discussed  "three  essential  parameters"  that 
define  a  locally  fitted  covariance  function,  namely,  the  variance,  the 
correlation  length,  and  the  curvature  parameter  of  the  covariance 
function,  after  the  function  is  restricted  to  the  sphere  where  the  given 
data  refer.  A  good  fit  to  the  variance  of  the  data  is  important  for 
yielding  realistic  error  estimates  from  the  collocation  procedure.  The 
correlation  length  carries  information  on  the  "smoothness"  of  the  field, 
i.e.,  on  the  spectral  distribution  of  the  total  variance.  The  curvature 
parameter,  measuring  the  curvature  of  the  covariance  function  at  V=0*, 
carries  information  about  the  very  local  behavior  of  the  field. 

A  discussion  on  the  determination  of  the  three  essential  parameters 
of  a  covariance  function  from  real  data  is  given  in  a  paper  by  Schwarz 
and  Lachapelle  (1980).  In  our  present  study,  where  the  data  were 
assumed  errorless  in  order  to  concentrate  on  modeling  errors,  a 
properly  fitted  variance  was  unimportant.  Also  in  our  study,  a  fit  to 
the  curvature  parameter  and  correlation  length  was  not  explicitly 
performed  unlike  in  the  above  paper.  Rather,  the  analytical  expression 
for  the  covariance  was  "tailored"  to  real  data  by  the  criteria  of  yielding 
the  least  RMS  error  of  actual  predicted  values  at  withheld  test  points  in 
the  local  area  (see  Table  8). 

For  the  purpose  of  covariance  function  modeling  it  is  useful  to  treat 
the  anomalous  gravity  field  as  being  generated  by  a  white  noise 
u«-layer  on  the  internal  sphere  of  radius  ro.  Actual  generations  of 
covariance  functions  from  white  noise  distributions  at  depth  is  discussed 
in  Heller  and  Jordan  (1979).  The  covariance  function  of  a  white  noise 
u»-layer  is  (Cruz  and  Laskowski,  1984,  sec.  6.2): 


cov(u*(P),  u*(J))  =  <To*  }  (2n+l)  P„(cos  Vpj)  .  (5.72) 

n=2 


where  <t%  is  the  variance  of  a  single  harmonic  uUnj.  On  the  other  hand, 
the  covariance  function  of  the  disturbing  potential  T*  on  the  ro -sphere 
is,  from  (5.70): 


cov(T»(P),  T*(J))  =  ro’  YL  (n-lp 

n=2 


Comparing  (5.72)  and  (5.73)  and  using  the  covariance  propagation 
concept,  then  we  can  conclude  that  T*  and  u*  are  spectrally  related  as 
follows: 
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!„♦  = 


n-1  ffo 


u„* 


(5.74) 


Equation  (5.74)  defines  the  u*  in  relation  to  the  disturbing  potential, 
and  is  analogous  to  equation  (5.23)  of  the  Dirac  approach.  For  example, 
if  d„=<ro  is  chosen  then  the  u*  is  seen  from  (5.74)  to  be  the  gravity 
anomaly  Ag«;  in  this  case  the  anomalous  gravity  field  can  be  interpreted 
as  being  generated  from  a  white  noise  gravity  anomaly  distribution  on 
the  internal  sphere.  From  the  discussions  of  Sjoberg  (1975,  1978),  one 
can  also  conclude  an  alternative  interpretation  of  the  u*  of  (5.74), 
namely,  that  the  (least  squares  collocation)  solution  essentially  minimizes 
the  variance  (sum  of  degree  variances)  of  the  quantity  u*  on  the 
Bjerhammar  sphere  of  radius  ro.  A  similar  interpretation  in  terms  of 
norm  minimization,  presented  in  a  different  way,  may  be  found  in 
Tscherning  (1972,  1983a). 

In  this  report  we  will  numerically  experiment  on  two  types  of 
choices  for  the  d„  in  (5.74).  These  choices  will  be  those  that  imply  the 
presence  on  the  Bjerhammar  sphere  of  a  white  noise  distribution  of  (1) 
gravity  anomaly  Ag«,  and  (2)  disturbing  potential  T*.  More  details  on 
these  two  systems  will  now  be  given. 

1.  the  white  noise  &g<-ay8teitt.  As  already  mentioned  the  choice  of 
dn  in  this  case  is: 


<i„(Ag*)  =  <7o. 


(5.75) 


The  covariance  function  (5.67)  then  becomes  exactly  in  the  form  of  the 
base  function  (5.18)  of  the  Dirac  approach,  except  for  the  replacement 
(5.68).  The  closed  covariance  expression  therefore  follows  easily  using 
(5.26)  and  the  replacement  (5.68): 


cov(T(P),  Agj)  =  Tp  s*  <ri  -t-l-SE-  s  cosT^  [s  +  3ln  ^  ^  cos'jH-E|  j  ^ 

(5.76) 


where 

r»» 

s  =  — 

rpFj 


(5.77) 


B*  =  1  -  2  s  cos  tC  s*  . 


•“a'/ 


.'•\v 


* 


(5.78) 


»T.V\' 
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Equation  (5.76)  checks  with  the  result  of  Sjoberg  (1975)  p.  106). 
Similarly)  we  can  use  the  replacement  (5.68)  into  (5.44)  -  (5.46)  to  obtain 
the  propagated  covariance  functions: 


cov(4g(P),  Agj)  =  <ri  ^  ! 


(5.79) 


cov(«,.(P).  Agj)  =  <Tg  s»  +  I  +  1  -  6E  - 

-  s  cos  V  (l3  +  6in  i-s_£2|3LLJj ]  .  (5.80) 


cov(tf^(P),  Agj) 

cos  a 

= 

s*  sin  i>  * 

cov(«\(P))  Agj) 

sin  a 

^  _  g  _  3  1 

-  s  cos'i'  -  E 

O  ,n  ^  *  COST^  +  E 

+  g  -  8  -  3  - 

E  sin* 

3  In  2 

(5.81) 


The  above  equations  all  agree  with  ibid.)  p.  106.  The  propagated 
covariance  functions  80*6  to  be  used  in  (5.66a)  for  computing  our 
functionals  of  interest. 

2.  the  white  noise  T*-system.  In  this  system  where  we  have  u*=T*) 
(5.74)  says  that  the  choice  of  d,,  must  be: 


d„(T*)  = 


(5.82) 


The  closed  covariance  expressions  of  interest  to  us  can  be  obtained  from 
Sjoberg  (ibid.)  p.  107): 


cov(T(P),  4gj)  =  <7?  -  5  +  s>]  ;  (5.83) 

cov(T(P),  T.)  =  <rg  ; 


(5.84) 


!.v»ji'  ■■)»  r,«  ’X  *-'  *  V  ».u  ".I,  -.w  »jr».i  ■.«  "5  "  J  *."  r.*  V  ."wr j  ’j.V-ny  rj  "j.  ^  t-jv.' 
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eov(A*(P,.  4gp=,|j^[i=ii|fii:-2eUz£ll*25^]  ,  ,53,, 


.  ri5{l-s»)^  12(2+s")(l-s»)  .  5+3s*l 

cov(«r(P).  A«j)  -  <rl  ^  [  -  ^-gV - ^ - |i -  +  — ; 


(5.86) 


i 


cov(fif(P),  Agj) 

cov(ax(P).  Agj) 


=  <to 


cos  a 


sin  a 


3  3^  sin  V'oj  r5(l-s*)* 
2  To*  E*  I  E» 


(5-s»)]- 

(5.87) 


The  white  noise  Tt-system  is  the  same  as  the  "attenuated  white  noise" 
covariance  model  of  Heller  and  Jordan  (1979)  up  to  a  scale  factor.  The 
said  covariance  model  was  recently  applied  by  Sunkel  (1983a)  in  his 
least  squares  collocation  modeling  of  the  geoid  in  Austria. 


This  completes  our  description  of  the  two  covariance  systemsi  the 
Ag*-  and  the  Tt-systems,  that  we  will  use  in  our  numerical  experiments 
on  the  least  squares  collocation  approach.  A  variety  of  other  covariance 
systems  may  be  defined,  by  simply  choosing  different  sets  of  d„-values 
provided  the  resulting  summations  converge  (see,  e.g.,  (5.69)).  It  turns 
out  that  covariance  systems  which  have  dp-values  of  the  same  "order  in 
n"  (see  second  paragraph  after  (5.39)  for  the  definition  of  "order  in  n") 
have  the  same  numericEil  characteristics.  This  is  the  reason  why,  for 
example,  models  of  degree  variances  (degree  variances  are  directly 
related  to  dp)  see  (5.69)  -  (5.71a))  have  been  distinguished  in  terms  of 
the  models’  variation  with  n  as  n  goes  to  infinity  (Tscherning,  1976,  p. 
2;  Moritz,  1980,  pp.  181-186).  Therefore,  the  two  covariance  systems 
that  we  have  chosen  for  experimentation  actually  typify  the  numerical 
behavior  of  two  groups  of  covariance  systems.  The  Agt-system  typifies 
systems  with  order  of  dp  equal  to  zero  (see  (5.75)),  and  the  T*-8y8tem 
typifies  systems  with  order  of  dp  equal  to  unity  (see  (5.82)). 

We  should  mention  that  the  well-known  Tscherning/Rapp  (1974) 
anomaly  degree  variance  model,  namely: 


/f: 


,  ^  A(n-1}  , 

"  (n-2)(n+24) 


(5.88) 


when  compared  with  (5.69)  yields  a  dp  with  order  equal  to  -1.  For  this 
model  the  only  effective  way  to  produce  the  very  short  correlation 
lengths  needed  for  local  computations  is  to  subtract  lower  harmonics 


3<n<N  from  the  covariance  functions  (ibid.,  p.  62).  This  removal  of 
lower  harmonics,  say  for  N=180,  is  expensive  in  terms  of  computer  time. 
For  this  reason  we  have  avoided  the  use  of  this  model.  This  is  the 
same  problem  that  caused  Sunkel  (1983a)  to  switch  to  the  use  of  Heller 
and  Jordan’s  (1979)  "attenuated  white  noise  model",  which  is  essentially 
the  T*-system  we  have  discussed  above.  For  the  Ag*-  and 
T*<-covariance  systems  that  we  will  use,  the  correlation  length  is 
adjusted  simply  by  changing  the  radius  ro  of  the  Bjerhammar  sphere  in 
the  closed  covariance  expressions. 

Finally,  we  should  discuss  the  computation  of  the  collocation 
coefficients  bj  (J=l,  2,  ...,  k)  appearing  in  (5.65).  As  in  the  Dirac 
approach  we  will  solve  the  coefficients  from  k  given  gravity  anomaly 
values  on  the  surface  of  the  earth.  We  first  form  the  k  x  k  system  of 
linear  equations  from  (5.66a): 


Ag,  =  bj  cov(Ag,,  Agj),  i 

j  =  l 


12  k 


(5.89) 


where  in  our  case  cov(Agj,  Agj)  will  be  either  (5.79)  or  (5.85).  Equation 
(5.89)  is  analogous  to  (5.58)  of  the  Dirac  approach.  For  a  direct 
comparison  between  the  Dirac  and  least  squares  collocation  approaches 
we  will  also  employ  the  Gauss-Seidel  iteration  to  solve  (5.89).  (Note  that 
Katsambalos  (1981)  used  the  direct  matrix  inversion  to  solve  (5.89),  but 
employed  the  Gauss-Seidel  iteration  to  solve  (5.58)).  The  employed 
Gauss-Seidel  algorithm  has  already  been  given  in  (5.59)  -  (5.64),  the 
only  change  being  the  replacement  of  (5.60)  by: 


A,  j  =  cov(Ag, ,  Agj) 


(5.90) 


This  completes  our  description  of  the  least  squares  collocation 
approaches  to  the  collocation  problem. 


5.4  Numerical  Investigations 


5.4.1  Tailoring  of  Covariance  Functions 


The  theory  of  least  squares  collocation  requires  that  the  covariance 
function  being  used  approximates  the  local  empirical  covariance  function 
in  the  area  of  computations.  This  "tailoring"  or  fitting  of  an  analytical 
covariance  expression  to  an  empirical  one  is  normally  performed  with 
respect  to  the  gravity  anomaly  auto-covariance  function.  Also,  it  is 
usually  sufficient  to  fit  just  the  variance  and  correlation  length, 
ignoring  the  very  irregular  curvature  parameter.  A  fit  to  the  variance 
of  local  data  does  not  present  any  problems,  as  this  can  be  done  by 


simply  scaling  the  entire  analytical  covariance  function  to  reflect  the 
desired  vEuriance  (value  at  ii-0').  A  fit  to  the  correlation  length,  on  the 
other  hand,  will  be  simple  if  it  can  be  done  by  Just  changing  the  radius 
of  the  Bjerhammar  sphere  occurring  in  the  covariance  expression.  This 
is  indeed  the  case  for  the  Ag*  and  T*  covariance  systems  that  we  have 
chosen  for  our  numerical  experiments. 

Figure  26  shows  the  covariance  functions  on  the  mean  earth  sphere 
R  for  the  Ag*  and  T*  systems,  as  the  depth  D  to  the  Bjerhammar  sphere 
takes  on  values  D=10,  20,  30,  50,  100,  500  km.  The  functions  were  scaled 
to  a  variance  of  1  mgal*.  Figure  26  resulted  from  (5.79)  for  the 
Agt-system  and  (5.85)  for  the  T»-sy8tem.  For  the  quantity  s  in  (5.77) 
we  set  rp=rj=R=6371  km,  and  ro=R-D.  The  graphs  reveal  that  to  a  good 
approximation  the  correlation  length  (  of  the  Ag*-system  is  related  to  D 
by: 


f(Ag*)  «  1.5  D  . 


(5.91) 


SlDilarly,  we  find  for  the  T*-syste*: 


((T*)  a  0.75  D  . 


(5.92) 


Equation  (5.91)  is  still  true  if  Ag*  is  replaced  by  other  white  noise 
systems  with  0(dn)=0  (see  discussion  above  (5.88).  In  fact,  we  found 
practically  the  same  covariance  graphs  when  using  systems  generated 
by  white  noise  distribution  on  the  internal  sphere  of  single  density 
layer,  potential  gradient,  radial  disturbance,  and  total  deflection  of  the 
vertical,  all  these  systems  having  0(d„)=0.  The  covariance  expressions 
for  these  systems  are  given  in  Sjoberg  (1975).  On  the  other  hand, 
(5.92)  remains  valid  also  for  the  covariance  system  generated  from  white 
noise  double  layer  distribution  on  the  Bjerhammar  sphere,  this  system 
having  0(dn)=l  like  the  T**-8yatem.  Again,  covariance  expressions  based 
on  a  white  noise  double  layer  distribution  may  be  found  in  Sjoberg 
(1975).  The  point  is  that  for  the  described  systems  of  0(d„)=0  and 
0(d„)  =  1  a  desired  correlation  length  (  may  be  easily  implemented  by 
specifying  the  proper  value  for  D.  Specifically,  very  short  (-values  on 
the  order  of  15  km,  needed  in  detailed  local  applications,  can  be  easily 
set. 


In  contrast,  the  Tscherning  and  Rapp  (1974)  model  entails  a  more 
difficult  tailoring  procedure  to  reach  very  short  correlation  lengths,  as 
Figure  27  shows.  The  only  effective  way  to  reach  (  values  on  the  order 
of  15  km  is  to  remove  lower  harmonics  3*n*180,  which  is  an  expensive 
procedure.  Changing  the  radius  of  the  Bjerhammar  sphere  helps  the 
tailoring  procedure,  but  only  in  a  very  limited  way  since  the  use  of 
small  D  (e.g.,  D=1  km)  is  still  not  enough  to  reach  (  values  on  the  order 
of  15  km. 
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Figure  26.  Gravity  Anomaly  Auto-Covariance  Functions  on  the  Mean 
Earth  Sphere  (R=6371  km),  Generated  from  White  Nuiae 
Distribution  of  Gravity  Anomaly  (top  diagram)  and  Disturbing 
Potential  (bottom  diagram)  on  an  Internal  Sphere  at  Depth 
D=10,  20,  30,  50,  100,  500  km.  Diagrams  reveal  that  for  the 
white  noise  gravity  anomaly  system  there  is  correlation 
length  (  ^  1.5  D,  and  for  the  white  noise  disturbing 


potential  system  (  -  0.75  D. 
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Figure  27.  Gravity  Anomaly  Auto-Covariance  Functions  on  the  Mean 
Earth  Sphere  (R=6371  km),  Tailored  from  the  Tscherning  and 
Rapp  (1974)  Covariance  Model  by  Using  Different  Depths  D 
to  the  Internal  Sphere  (D=l,  2,  3,  5,  10,  15  km)  and 
Subtracting  Lower  Harmonics  3<n<N  (left  to  right:  N=20,  36, 
180).  Very  short  correlation  lengths  on  the  order  of  15  km 
can  only  be  obtained  by  subtracting  the  lower  180 
harmonics. 
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5.4.2  Real  Data  Predictions  of  Gravity  Anomalies  in  Canada 

Sjober^  (1978,  Tables  8.1  and  8.2)  gives  surface  gravity  anomaly 
data  with  station  elevations,  that  can  be  used  for  intercomparison  of 
gravity  anomaly  prediction  techniques.  The  points,  located  in  Manitoba, 
Canada,  are  shown  in  Figure  28.  Eighty-seven  points  (shown  as  circles) 
are  considered  as  "observation  points"  and  50  points  (shown  as  crosses) 
are  considered  as  "prediction  points"  but  are  also  known.  An  important 
characteristics  of  the  point  distribution  is  that  the  observation  points 
are  spaced  almost  regularly  in  a  0t5x0‘.5  grid,  while  the  prediction  points 
are  generally  near  the  centers  of  the  data  gaps.  Since  the  actual 
gravity  anomaly  function  can  easily  have  features  beyond  the  maximum 
0’.5x0‘.5  resolution  implied  by  the  data  points,  one  should  not  expect  a 
highly  accurate  prediction  of  the  withheld  test  points. 

The  non-uniqueness  of  the  collocation  solution  in  between  the  data 
points  is  illustrated  in  Figure  29.  The  two  collocation  solutions  shown 
are  for  the  2*x3*  area  marked  off  in  Figure  28.  The  solutions  were 
generated  using  the  Dirac  approach  with  gravity  anomaly  impulses  on 
the  Bjerhammar  sphere  of  radius  =  (6371  km-D).  The  depth  D=30  km 
was  used  for  the  solution  shown  on  the  left,  while  D=75  km  was  used  for 
the  solution  on  the  right.  In  accordance  with  the  collocation  condition 
(5.9)  the  two  solutions  agree  closely  with  the  data  values  and  with  each 
other  at  the  approximately  0:5xD’.5  grid  locations  where  the  data  are 
given.  The  agreement  at  the  data  points  is  not  exact  in  this  case, 
because  of  the  use  of  an  iterative  linear  equation  solver  (see  (5.59) 
which  was  terminated  aifter  5  iterations.  Although  the  solutions 
approximately  agree  at  the  data  locations,  they  may  have  very  different 
behavior  in  between  the  data  points.  As  expected  the  D=75  km  solution 
appears  smoother  (larger  correlation  length)  than  the  D:30  km  solution. 
It  is  important  to  state  that  both  solutions  were  contoured  from 
predicted  point  values  in  a  0t25x0’.25  grid,  not  just  a  Ot5xO'.5  grid.  This 
was  to  make  sure  that  the  surface  shown  in  between  the  original  0:5x0*.5 
data  points  was  the  actual  predicted  surface,  and  not  the  surface  which 
would  have  been  artificially  produced  by  the  contouring  algorithm  if 
contours  were  generated  from  0'.5x0’.5  point  values. 

Which  then  is  the  optimal  Bjerhammar  sphere  radius  to  use  under  a 
given  base  function  system  in  collocation?  Different  techniques  have 
been  devised  in  the  literature  to  answer  this  question.  The  reader 
interested  in  details  should  consult  Needham  (1970),  Blaha  (1983),  and 
Sunkel  (1981b,  1983b)  for  the  case  of  the  point  mass  system,  which  is 
sdso  applicable  to  any  Dirac  system.  For  an  optimal  internal  sphere 
radius  in  the  least  squares  collocation  approach  and  for  the  more 
general  minimum  norm  collocation  approach,  see  Sjoberg  (1975)  and 
Lelgemann  (1981).  For  the  least  squares  collocation  approach,  there  is  a 
built-in  theory  to  determine  an  optimal  radius,  namely,  to  use  a  radius 
that  will  cause  the  covariance  function  to  fit  the  correlation  length  of 
local  data  in  the  area  of  computation.  However,  the  most  direct  method 
to  determine  an  optimal  radius  from  the  data  themselves  is  to  examine 
the  prediction  errors  at  withheld  test  points  in  the  area,  for  various 
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Point  Location  of  Surface  Gravity  Anomalies  in  Manitoba, 
Canada  (plotted  from  Sjoberg,  1978,  Tables  8.1  and  8.2). 
(Circles:  observation  points;  crosses:  withheld  prediction 

points). 
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•e  29.  Two  Gravity  Anomaly  Solutions  (C.l.=5  mgals)  for  the  \:cix 
marked  off  in  Figure  28.  Predicted  Using  the  Dirac 
Approach  with  Gravity  Anomaly  Impulses  on  the  Internal 
Sphere  of  Radius  =  (6371  km-D).  Left:  D=30  km;  Right: 

D=75  km;  No.  of  iterations=5. 
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choices  of  internal  sphere  radius.  This  last  approach  was  specifically 
applied  by  Sjoberg  (1975,  1978)  in  order  to  determine  an  optimal  radius 
for  a  given  base  function.  It  is  here  that  Sjoberg  concluded  that  as 
long  as  an  optimal  radius  is  used,  various  base  function  systems  yield 
the  same  quality  of  results. 

We  decided  to  repeat  the  gravity  anomaly  predictions  of  Sjoberg 
(1978,  Table  8.3)  since  conflicting  results  were  presented  by  Katsambalos 
(1981,  Table  9.3).  The  data  used  were  the  observation  points  and 
withheld  teat  pxDints  shown  in  Figure  28.  The  methods  compared  were 
the  least  squares  collocation  approach  with  white  noise  gravity  anomalies 
on  the  internal  sphere,  which  we  denote  by  l.a.c.  (4g*),  versus  the  Dirac 
approach  with  gravity  anomaly  impulses  on  the  internal  sphere,  which 
we  denote  by  Dirac  (ig*).  We  additionally  used  the  other  two  Dirac 
systems  that  we  discussed  earlier,  namely,  the  point  mass  and  point 
dipole  Dirac  systems,  which  we  denote  by  Dirac  (>i*)  and  Dirac(7i*), 
respectively.  The  number  of  iterations  to  satisfy  the  criteria  (5.64),  and 
RMS  prediction  error  at  the  50  test  points,  are  shown  in  Table  8  for 
various  depths  D  to  the  internal  sphere. 

Our  results  are  in  complete  agreement  with  Sjoberg’s  Table  8,3,  for 
the  l.s.c.  (4g*)  and  Dirac  (4g»)  solutions.  At  the  optimal  depths  of  15 
km  for  the  l.s.c.  30  km  for  the  Dirac  (4g*),  30  km  for  the  Dirac 

(p*),  and  50  km  for  the  Dirac  (m*),  the  various  systems  yielded  about 
the  same  RMS  error  of  10.2  mgal,  and  required  about  the  same  number 
of  iterations  of  5  for  convergence  under  criteria  (5.64).  Since  the  RMS 
anomaly  for  the  50  prediction  points  is  13.53  mgals,  we  have  the 
following  normalized  measure  of  "unrecovered"  information  (see 
Tscherning,  1981,  sec.  4): 


RMS  Variation  (Observed  -  Computed) 
RMS  Variation  of  Observations 


(5.93) 


10.2 

13.5 


76%  , 


which  is  within  expectation,  considering  the  large  data  spacing  of  0'.5x0'.5 
(see  ibid..  Figure  1). 

Table  8  shows  that  the  Dirac  (Ag*)  and  Dirac  (/i*)  systems  produced 
practically  identical  results,  a  consequence  of  the  fact  that  both  systems 
have  0(kn)=0  (see  discussion  above  (5.40)).  The  Dirac  (7i*)  system  which 
has  0(kn)  =  i  behaved  differently,  requiring  a  deeper  Bjerhammar  sphere 
than  the  other  two  Dirac  systems.  As  observed  by  Sjoberg  (1978)  and 
proven  by  him  theoretically,  the  system  of  equations  arising  from  the 
Dirac  formulation  is  inherently  more  stable  than  that  arising  from  the 
least  squares  collocation  approach,  for  a  given  radius  of  the  internal 
sphere.  This  difference  in  stability  is  reflected  in  Table  8  as 
differences  in  required  number  of  iterations.  .Also,  the  l.s.c.  solution  is 
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Table  8.  "Number  of  Iterations/RMS  Error  (mgal)"  Under  Various 
Discrete  Approaches,  When  Predicting  50  Withheld  Gravity 
Anomalies  from  87  Known  Gravity  Anomalies  in  Manitoba, 
Canada. 


Depth  to 
Internal 
Sphere()an) 

l.s.c.  (Ag*) 

Dirac  (Ag*) 

Dirac  (p*) 

— 

0 

■m 

1/13.5 

5 

1/13.6 

10 

2/12.0 

2/13.6 

15t 

5/10. 2t 

2/11.0 

2/13.1 

20 

11/10.4 

3/10.4 

3/10.4 

2/12.5 

25 

25/10.6 

3/10.2 

3/10.2 

2/11.9 

30t 

30/1Q.9 

4/10. 2t 

4/10. 2t 

2/11.4 

35 

30/11.1 

5/10.2 

5/10.2 

3/10.9 

40 

30/11.3 

7/10.3 

7/10.3 

3/10.6 

45 

30/11.5 

9/10.5 

9/10.5 

4/10.5 

50t 

30/11.7 

12/10.6 

12/10.6 

5/10. 5t 

55 

30/11.8 

16/10.7 

16/10.7 

6/10.5 

60 

8/10.6 

65 

10/10.7 

70 

14/10.8 

75 

18/10.9 

80 

24/11.0 

85 

30/11.2 

90 

30/11.3 

t  optimal  results.  Iterations  were  terminated  when  RMS  Residual  <  0.25 
mgal,  or  Maximum  Residual  <  0.5  mgal,  or  No.  of  Iterations  >  30. 


much  more  sensitive  than  the  Dirac  solution,  to  changes  in  the  radius  of 
the  internal  sphere.  This  is  evidenced  by  the  narrow  range  of  useful 
depths  (e.g.,  10  km  to  25  km)  that  can  be  used  for  the  l.s.c.  method,  as 
compared  with  the  broad  range  of  useful  depths  for  the  Dirac  systems 
(20  km  to  55  km  for  the  Dirac  (Ag*)  and  Dirac  30  km  to  65  km  for 

the  Dirac  (p*)). 


5.4.3  Comparison  of  Discrete  and  Continuous  Models  in  the  New  Mexico 


Test  Area 


In  this  section  we  compare  the  various  models  discussed  earlier,  in 
terms  of  their  ability  to  model  the  high  frequency  portion  of  the  radial. 


North-South,  and  East-West  components  of  the  disturbance  vector  in 
space.  The  problem  is  a^ain  the  computation  of  the  external  field  from 
gravity  anomaly  data  on  the  earth’s  (non-spherical)  topography.  We 
concentrate  on  the  high  frequency  part  of  the  total  field,  because  this 
is  the  part  affected  most  by  the  topography  problem.  We  will 
intercompare  the  discrete  solutions  (solutions  to  the  Bjerharamar 
problem)  discussed  in  this  chapter,  as  well  as  the  classical  integral 
solution  used  in  Chapter  4.  Computation  points  will  be  along  the  five 
vertical  test  lines  in  a  l*xl*  area  in  New  Mexico,  already  described  in 
Chapter  4  (Figure  25).  Altitudes  of  computation  will  again  be  at  5,  10, 
30,  100,  and  500  km,  giving  25  test  points  in  all  (5  points  per  line,  times 
5  test  lines). 

To  serve  as  the  "true"  field  we  took  the  residual  topographic  model 
(RTM),  generated  from  the  gravitational  influence  of  the  positive  and 
negative  topographic  masses  lying  between  the  actual  topography  and  a 
reference  topography  to  spherical  harmonic  degree  180.  Specifically,  we 
used  the  field  generated  in  Test  7  of  Chapter  4.  The  field  was 
generated  from  topographic  data  in  a  2*x2*  area  centered  at  the 
computation  foot  points.  The  integration  of  gravitational  influences  was 
performed  using  Forsberg's  (1984)  prism  integration  program  along  with 
30"x30"  point  elevations,  5'x5'  mean  elevations,  reference  topographic 
elevations  to  degree  180,  and  an  assumed  density  of  2.67  g/cm’.  Values 
were  generated  for  5'x5'  surface  gravity  anomalies  in  the  2*x2*  area 
and  for  the  disturbance  components  at  the  25  test  points  in  space. 
These  values  were  then  considered  as  the  self-consistent  or  "true"  set 
of  values  for  testing  different  techniques  for  modeling  the  spatial 
disturbance  vector  field  from  surface  gravity  anomaly  data.  The  three 
disturbance  components  at  each  of  the  25  test  points  in  space  are  shown 
in  Table  9.  The  overall  RMS  of  the  values  in  Table  9  is  12.36  mgals. 

As  a  first  step  we  performed  computations  analogous  to  Table  8  to 
determine  for  each  base  function  system  an  optimal  internal  sphere 
radius  to  use  with  the  given  surface  gravity  anomaly  data.  This  time, 
optimal  radii  were  determined  not  by  minimizing  the  RMS  anomaly 
prediction  error  at  withheld  surface  points,  but  by  minimizing  the  RMS 
difference  between  computed  and  "true”  (Table  9)  disturbance 
components  in  space.  The  results  are  shown  in  Table  10,  giving  the 
number  of  iterations  performed  (see  (5.59))  and  the  RMS  disturbance 
error  (computed  -  "true"  value)  for  various  depths  to  the  Bjerhammar 
sphere  and  various  base  function  systems.  To  repeat,  for  visualization 
one  can  use  the  interpretations  that:  the  Dirac  (4g*)  solves  for  gravity 
anomaly  impulses  on  the  internal  sphere;  the  Dirac  (p*)  solves  for  point 
dipoles  on  the  internal  sphere;  the  l.s.c.  (4g*)  uses  the  attenuated  white 
noise  gravity  anomaly  covariance  function;  and  the  l.s.c.  (T*)  system 
uses  the  attenuated  white  noise  potential  covariance  function.  The 
number  of  iterations  allowed  was  reduced  from  the  30  used  in  Table  8  to 
just  10,  as  the  results  did  not  significantly  improve  with  more 
iterations.  The  point  mass  system  Dirac  (;j*)  was  confirmed  to  yield 
practically  identical  results  as  the  Dirac  (Ag*)  system,  and  is  now 
omitted  from  our  numerical  presentations.  Based  on  the  results  of 


Table  10  we  selected  the  optimal  depth  of  10  km  for  the  Dirac  (4g*),  15 
km  for  the  Dirac  {]!*),  3.5  km  for  the  l.s.c.  (Ag*),  and  15  km  for  the 
l.a.c.  (T*).  Note  that  for  the  Dirac  (Ag*)  the  result  for  depth  20  km  is 
even  slightly  better  than  that  for  depth  10  km,  but  a  shallower  depth  is 
preferable  because  in  principle  the  system  will  be  more  stable.  The 
results  corresponding  to  these  optimal  depths  were  then  used  in  any 
further  comparisons  of  the  different  methods.  Note  from  Table  10  that 
the  Dirac  systems  have  faster  convergence  (fewer  number  of  iterations) 
than  the  least  squares  collocation  systems.  Also,  it  is  indicated  that  the 
Dirac  systems  have  better  agreement  with  the  residual  topographic  model 
(RTM)  than  do  the  l.s.c.  systems. 


Table  10.  "Number  of  Iterations/RMS  Error  (mgals)"  Under  Various 
Discrete  Approaches,  When  Modeling  the  Disturbance 
Components  of  the  Residual  Topographic  Field,  from  5'x5' 
Surface  Gravity  Anomalies  in  the  2*x2*  New  Mexico  Test  Area. 


Depth  to 
Internal 
Sphere  (km) 

Dirac  (Ag*) 

Dirac  (/i*) 

l.s.c.  (Ag*) 

l.s.c.  (?*' 

1 

3/6.01 

2/11.88 

10/2.39 

10/12. 17 

3.5 

5/1.79 

3/8.50 

10/1. 55t 

10/6.28 

5 

7/1.23 

3/5.21 

10/1.96 

10  '4.12 

10 

10/1. 14t 

5/1.08 

10/2.57 

10/2.26 

15 

10/1.23 

10/1. 05t 

10/3.17 

10/1. 82t 

20 

10/0.98 

10/1.11 

10/2.96 

10/2.67 

t  optimal  results.  Iterations  were  terminated  when  RMS  Residual  <  0.25 
mgal,  or  Maximum  Residual  <  0.5  mgal,  or  No.  of  Iterations  >  10. 

Next,  we  now  attempt  to  provide  a  specific  feeling  for  the  numerical 
differences  among  the  various  discrete  models,  and  the  classical  integral 
model  used  in  Chapter  4.  As  presented  in  Chapter  4,  the  integral  model 
can  account  for  the  mean  topography  through  the  quantity  hm  in  the 
formulas,  where  h„  is  a  mean  topographic  elevation  in  the  area  of 
computation.  Obviously,  the  neglect  of  topographic  variations  around 
the  mean  is  a  modeling  error  being  committed  in  this  approach.  Our 
present  comparisons  will  give  a  feeling  for  the  magnitude  of  this  error 
(in  fact,  this  error  was  also  examined  in  Chapter  4).  The  interest  in 
using  the  integral  model  lies  in  its  relative  economy,  not  requiring  any 
solution  of  a  linear  equation  system  in  contrast  to  the  discrete  models. 
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We  compared  the  various  models  by  pairs,  resulting  in  a  total  of 
nine  pairings.  In  every  pairing,  each  of  the  two  models  involved  had 
its  own  set  of  75  disturbance  values  analogous  to  Table  9.  We 
differenced  the  two  tables  of  disturbances  involved  in  each  pairing, 
producing  nine  tables  containing  75  disturbance  differences  each.  Based 
on  looking  at  the  nine  tables  of  differences  we  decided  that  for  the 
analysis  and  presentation  of  results  it  would  be  sufficient  to  represent 
each  table  by: 


•  the  RMS  value  of  the  75  differences  contained  in  the  table;  and 

•  the  test  line  or  lines  (out  of  A,  B,  C,  D,  E;  see  Chapter  4)  which 

show  maximal  differences  for  the  6^,  5^,  fix  disturbance 

components. 


The  results  for  the  nine  pairings  are  shown  in  Table  11. 

For  the  integral  model  involved  in  pairings  (1)  and  (9),  a  mean 
topographic  elevation  of  h„=1.5  km  was  used  as  in  Chapter  4,  When  the 
value  of  h„=0  km  was  tested,  thereby  assuming  the  given  gravity 
anomalies  to  refer  to  the  geoid,  the  RMS  difference  with  the  RTM 
(pairing  (1))  increased  significantly  from  1.75  mgal  to  2.26  mgal,  and  the 
RMS  difference  with  the  Dirac  (4g*)  (pairing  (9))  increased  from  0.98 
mgal  to  1.65  mgal.  Therefore,  we  specifically  conclude  that  there  is  a 
significant  gain  in  accuracy  when  the  mean  topography  is  accounted  for 
through  h„  in  the  classical  integral  formulas. 

Pairings  (1)  -  (5)  show  the  comparison  of  the  integral  model  plus 
the  four  "discrete  approach"  models  against  the  residual  topographic 
model  (RTM).  Again,  we  assume  the  RTM  to  be  the  "true"  field  and  call 
the  differences  shown  as  "errors". 

Comparing  (1)  with  either  (2)  or  (3),  one  can  see  a  significant 
although  not  drastic  increase  in  accuracy  with  the  use  of  the  Dirac 
model  over  the  integral  model.  Improvement  is  especially  evident  at 
altitudes  10  km  and  30  km.  Above  30  km  the  choice  of  model  is  not 
critical,  as  high  frequency  effects  are  greatly  attenuated  at  such 
altitudes.  The  residual  errors  of  the  Dirac  models  especially  at  low 
altitudes  can  be  attributed  to  the  use  of  5'x5'  gravity  anomalies  and 
elevations,  which  automatically  limits  the  resolution  of  field 
representation.  The  RTM  inherently  has  more  detailed  resolution,  with 
its  use  of  up  to  1  km  x  1  km  topographic  height  data  (which  directly 
converts  to  1  km  x  1  km  resolution  of  gravity  data  because  of  correltion 
between  elevation  and  gravity  anomalies).  Therefore,  we  conclude  a 
significant  gain  in  accuracy  with  the  use  of  a  Dirac  model  over  the 
integral  model.  But  at  the  same  time,  the  Dirac  model  (or  any  of  the 
other  collocation  models)  cannot  represent  the  very  high  resolution 
attainable  with  the  RTM  at  very  low  altitudes.  The  RTM  should  always 
be  used  as  one  of  the  complementary  models  at  low  altitudes,  whenever 


detailed  height  information  is  available. 

Examining  pairings  (4)  and  (5),  no  evidence  can  be  found  of  any 
improvement  with  the  use  of  the  least  squares  collocation  (l.s.c.)  models 
over  the  Dirac  or  the  integral  models.  The  problem  encountered  during 
computations  with  the  l.s.c.  systems  was  the  instability  of  linear 
equation  systems.  The  use  of  more  widely  spaced  gravity  anomaly  data 
(say,  lO'xlO')  would  have  made  the  systems  more  stable,  but  this  would 
also  have  automatically  reduced  the  achievable  resolution  of  the  field. 
We  conclude  that  for  very  detailed  field  representations,  the  stability  of 
the  arising  linear  equation  system  plays  a  significant  role  in  the 
discrete  approaches,  and  in  such  cases  the  Dirac  systems  are  preferable 
to  the  least  squares  collocation  systems. 

Examining  pairing  (6),  one  can  see  that  with  the  use  of  their 
respective  optimal  radii  to  the  internal  sphere,  the  Dirac  systems  agree 
with  each  other  in  the  submilligal  level  even  at  the  low  altitude  of  5  km. 
Therefore,  from  this  test  we  do  not  see  any  preference  to  the  point 
gravity  anomaly,  or  the  point  mass,  or  the  point  dipole  systems  with 
respect  to  each  other. 

On  the  other  hand,  pairing  (7)  shows  that  the  l.s.c.  system  is  much 
more  sensitive  to  the  choice  of  either  the  T*-  or  the  4g*-syatem  (This  is 
equivalent  to  the  choice  of  norm  to  be  minimized;  see  remarks  after 
(5.74)).  Also,  the  great  sensitivity  of  the  l.s.c.  systems  to  the  depth  of 
the  Bjerhammar  sphere  (See  Tables  10  and  8)  contributes  to  the  dif¬ 
ferences  found  in  pairing  (7).  Also,  the  results  for  the  l.s.c.  systems 
refer  to  linear  equation  solutions  that  have  not  satisfactorily  converged 
'see  (5.64)  for  convergence  criteria)  in  the  course  of  the  iterations. 
Therefore,  convergence  problems  also  play  a  role  in  the  differences 
found  in  pairing  (7).  We  should  note  that  direct  matrix  inversions  were 
not  tried  in  these  types  of  tests,  but  they  should  have  been. 

Pairing  (8)  compares  the  Dirac  (4g*)  with  the  l.s.c.  (dg*).  Again, 
the  significant  differences  can  be  attributed  to  the  stability  problems 
mentioned  above  for  the  l.s.c.  system. 

Finally,  pairing  (9)  shows  how  the  integral  model  performs  in  com¬ 
parison  with  the  Dirac  model.  As  one  can  see,  the  differences  between 
the  two  models  are  not  very  large,  although  they  are  significant.  The 
differences  indicate  how  much  additional  information  about  the 
disturbance  field  is  gained  by  accounting  for  the  full  variations  of  the 
topography  (as  in  the  Dirac  system),  not  just  the  mean  topography  (as 
in  the  integral  system).  Above  30  km  the  integral  model  is  sufficient, 
but  at  very  low  altitudes  the  Dirac  system  should  be  used  for  accurate 
determinations. 


Table  11.  RMS  and  Maximal  Differences  Between 
Disturbance  Vector  Components  (mgals). 


Various  Models 


RMS 

Difference 


1.75 

[Integral 
-  RTM] 


(2)  1.14 

[Dirac  (Ag*) 
-  RIM] 


(2)  1.05 
[Dirac  (JI*) 
-  RIM] 


1.55 

[l.s.c.  (Ag*) 
-  RTM] 


1.82 

[l.s.c.  (!♦) 
-  RIM] 


0.49 

[Dirac  (Ag*) 

-  Dirac  (Ji*)  ] 


1.38 

[l.s.c.  (Ag*) 
-  l.s.c.  (T*) 


(8)  , 


D  6r 

A 

D  Ax 


1.30  A  6^ 

[Dirac  (Ag*)  C  Ay 
-  l.s.c.  (Ag*)]  D  Ax 


0.98 

[ Integral 
-  Dirac  ( Ag*) 1 


Test  Lines  with  Maximal  Differences 


-4.07 

3.10 

-3.96 


-1.85 

-2.20 
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-3.38 

2.67 

-2.26 


-1.48 
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6.  SUMMARY,  CONCLUSIONS,  RECOMMENDATIONS 


In  this  report  we  have  discussed  the  modeling  of  the  external 
gravity  disturbance  vector  of  the  earth,  from  surface  gravity  anomaly 
data.  This  work  is  a  continuation  of  the  work  of  Katsambalos  (1981) 
with  important  differences  and  expansion  of  treatment  as  noted  in 
Section  1.1.  In  principle  the  disturbance  vector  is  modeled  as 

"S  -  grad  T,  with  the  disturbing  potential  T  being  in  turn  related  to  the 
gravity  anomaly  Ag  through  (1.9)  or  its  spherical  approximation  (1.12). 
Under  the  assumption  that  T  is  harmonic  outside  the  earth's  attracting 
masses,  the  knowledge  for  Ag  at  every  point  on  the  earth's  surface  is 
sufficient  to  define  T  in  all  of  apace.  More  details  on  preliminary 
concepts  are  given  in  Section  1.2,  discussing  the  precise  definition  of 
Ag,  the  relation  between  Ag  and  T  at  a  point,  and  the  realization  of  Ag 
from  observational  data  on  the  earth's  surface.  In  Section  1.3  a 
discussion  is  given  outlining  the  scope  of  the  entire  study. 

Chapter  2  discusses  the  re-parameterization  of  the  long-wavelength 
components  of  the  anomalous  gravity  field,  from  the  given  surface 
gravity  anomalies  to  a  set  of  spherical  harmonic  coefficients  of  the 
earth's  disturbing  potential.  The  focus  is  on  the  analytical  continuation 
of  the  surface  gravity  anomaly  data  to  values  that  refer  to  a  spherical 
boundary  surface.  The  Ag-values  on  a  sphere  are  the  values  directly 
useable  in  the  determination  of  spherical  harmonic  coefficients.  It  is 
recommended  that  the  analytical  continuation  be  carried  out  in  two 
steps:  (a)  analytical  continuation  of  surface  values  to  the  ellipsoid; 
followed  by  (b)  analytical  continuation  of  values  from  the  ellipsoid  to  a 
sphere,  moat  appropriately  the  equatorial  sphere. 

Step  (a)  above  can  be  implemented  by  Taylor  series  expansion  using 
the  vertical  gradients  of  the  field.  Operationally,  such  gradients  can  be 
obtained  from  an  existing  spherical  harmonic  expansion,  in  application  of 
the  ideas  of  Rapp  (1984).  Step  (b)  can  also  be  implemented  by  Taylor 
series,  but  because  the  deviation  between  ellipsoid  and  sphere  is  a 
simple  function  of  latitude  it  is  possible  in  this  case  to  analytically 
expand  the  terms  of  the  Taylor  series  into  spherical  harmonics. 
Specifically,  Taylor  terms  involving  up  to  second  order  gradients  have 
been  expanded  in  Chapter  2.  Our  main  contribution  via  Chapter  2  is  in 
showing  that  such  spherical  harmonic  expansions  of  Taylor  series  terms 
lead  to  the  ellipsoidal  correction  terms  of  Pellinen  (1982).  As  a  result,  it 
is  shown  that  the  application  of  Pellinen’s  ellipsoidal  corrections  is  an 
analytical  implementation  of  step  (b)  above.  This  is  in  contrast  to  the 
numerical  implementation  of  step  (a).  It  is  also  shown  that  the  use  of 
the  ellipsoidal  corrections  as  given  in  Pellinen  (ibid.)  implies  the  use  of 
a  relation  between  Ag  and  T  that  neglects  only  terms  of  0(e*),  i.e.,  the 
relation  (1.13). 
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Tables  1  and  2  of  Chapter  2  indicate  the  number  of  Taylor  series 
terms  that  need  to  be  considered  in  the  analytical  continuation  of  values 
from  ellipsoid  to  sphere.  The  higher  the  maximum  degree 
spherical  harmonic  series  development  the  more  Taylor  terms  need  to  be 
used.  For  example,  with  N„a„  =  300  it  is  indicated  that  terms  up  to  the 
second-order  gradient  need  to  be  considered,  leaving  (gravity  anomaly, 
height  anomaly)  errors  of  about  (0.8  mgal,  2  cm)  in  the  polar  areas. 
The  corresponding  effects  of  the  analytical  continuation  of  surface 
values  to  values  on  the  ellipsoid  are  an  order  of  magnitude  smaller, 
because  of  the  shorter  vertical  distances  involved  in  the  analytical 

continuation. 

For  the  complete  modeling  of  the  disturbance  vector  signal  in  its 
entire  frequency  range,  locally  valid  models  must  be  used  to  complement 
the  globally  valid  spherical  harmonic  model.  Before  studying  such  local 
models  in  Chapters  4  and  5,  however,  Chapter  3  first  gives  a 

familiarization  study  of  the  signal  characteristics  to  aid  the  design  of 
models  and  numerical  experiments.  Information  is  presented  on  the 
spectral  energy  distribution  (Section  3.2)  and  data  response 
characteristics  (Section  3.3)  of  the  disturbance  vector  signal  as  a 

function  of  altitude  in  space.  For  example,  according  to  Figures  4  and  5 

a  representation  to  n  =  360  (0*.5x0*.5  resolution)  will  resolve  at  least  90% 
of  the  RMS  signal  value  for  altitudes  above  30  km;  for  an  altitude  of  (5, 
10,  20)  km  such  representation  to  n  =  360  leaves  about  (16,  10,  5)  mgals 
unresolved. 

At  the  low  altitude  of  5  km,  for  example,  a  representation  to  n  = 
2500  (5'x5'  resolution)  is  needed  to  resolve  99%  of  the  RMS  radial  or 
horizontal  disturbance  signal.  On  a  point  by  point  computation  of  signal 
from  gravity  anomaly  data,  however,  it  is  not  necessary  to  have  data  of 
such  uniform  resolution  all  over  the  earth.  Rather,  data  of  less  and 
less  resolution  (and  accuracy)  can  be  used  farther  and  farther  away 
from  the  computation  point.  Such  data  response  characteristics  of  the 
signal  are  studied  in  Section  3.3  using  the  tools  of  truncation  theory. 
Under  the  formulation  of  truncation  theory  we  obtain  the  so-called 
truncation  error,  i.e.,  the  amount  of  error  propagated  into  the  computed 
signal  due  to  the  use  of  a  particular  set  of  data  in  the  "remote  zone". 
There  are  four  factors  affecting  the  truncation  error:  the  data 

accuracy,  the  data  resolution,  the  truncation  method  used  (i.e.,  the 
kernel  modification  used),  and  the  distance  of  the  data  from  the 
computation  point  (i.e.,  the  cap  size).  The  influences  of  these  factors 
are  shown  in  Figures  7  to  14,  for  the  radial  and  horizontal  disturbance 
components  at  a  high  altitude  of  100  km  and  a  low  altitude  of  5  km. 

As  an  example  of  the  use  of  Figures  7  to  14,  consider  first  Figure  7 
showing  the  case  of  radial  disturbance  computations  at  altitude  100  km. 
For  the  assumed  data  error  model  (see  (3.64)),  the  data  resolution  to 
Sr,f  =  20  (®  9*x9*  resolution),  the  unmodified  truncation  method,  and 
the  cap  size  Vo  =  13*,  we  find  from  Figure  7  that  the  propagated  error 
to  the  signal  is  about  0.2  ragal.  If  we  change  to  the  use  of  errorless 
data  (circled  plot  in  Figure  7)  the  required  cap  size  to  maintain  a  0.2 


mgal  error  reduces  from  Vo  =  13*  to  Vo  =  9*-  moreover,  we  increase 
the  data  resolution  to  Nref  =  36  (*  5*x5*  resolution),  then  the  required 
cap  size  for  a  0.2  mgal  error  reduces  further  from  Vo  =  9*  to  Vo  =  5*. 
Now  looking  at  Figure  8,  we  see  that  the  originally  required  cap  size  of 
Vo  =  13*  can  also  be  drastically  reduced  (to  Vo  =  6*)  if  we  change  from 
the  use  of  the  unmodified  truncation  to  the  use  of  the  improved 
Molodensky  truncation. 

General  conclusions  can  be  stated  based  on  Figures  7  to  14  as 
follows.  The  horizontal  disturbance  is  much  more  sensitive  to  remote 
zone  data  than  the  radial  disturbance,  as  evidenced  by  the  larger  data 
caps  (more  than  2  times  larger  than  the  radial  case)  needed  to  keep 
truncation  error  below,  e.g.,  a  0.2  mgal  level.  The  improved  Molodensky 
truncation  method  gives  better  results  than  the  the  unmodified 
truncation  for  relatively  large  cap  sizes  only;  anticipating  the  use  of 
remote  zone  data  to  resolution  n^ef  ^  180,  then  we  can  say  that  for  cap 
sizes  and  accuracy  levels  of  interest  the  improved  Molodensky  method 
offers  no  significant  gain  over  the  unmodified  method  for  radial  or 

horizontal  disturbance  computations  in  space. 

Chapter  4  starts  the  study  of  local  models  that  can  be  used  to 

complement  the  globally  valid  spherical  harmonic  model  studied  in 
Chapter  2.  In  the  global  discussions  of  Chapter  2,  the  problem  of  the 
shape  of  the  boundary  surface  lies  in  the  deviations  of  the  earth's 
surface  from  the  equatorial  sphere;  these  deviations  can  amount  to  20 
km  at  the  poles.  In  contrast,  for  the  local  models  such  as  those  of 
Chapter  4,  the  problem  with  the  shape  of  the  boundary  surface  consists 
only  in  the  deviations  of  the  earth's  surface  from  an  equipotential 
surface.  The  equipotential  surface  is  then  considered  as  a  sphere,  after 
a  formal  neglection  of  the  earth’s  flattening  under  the  spherical 

approximation  (see  Figure  2).  Therefore,  whereas  in  the  spherical 
harmonic  model  both  the  topography  and  ellipticity  are  important  issues 
in  the  modeling,  in  the  case  of  the  local  models  of  Chapters  4  and  5 
only  the  topography  needs  to  be  of  special  concern. 

In  Section  4.1  a  degree  variance  analysis  is  first  performed  to 
quantify  the  role  of  the  topography  in  the  modeling  problem.  The 

observed  degree  variances  of  the  earth’s  anomalous  gravity  field  are 
expressed  as  the  combined  effect  of  the  gravitational  influence  of 
condensed  topographic  masses  and  their  Airy-Heiskanen  isostatic 
compensation  at  depth  D.  It  is  then  numerically  shown  the  considerable 
smoothing  of  the  field  that  results  from  the  removal  of  shallow 
topographic  masses  associated  with  topographic  heights  beyond  harmonic 
degree  N  (see  Figures  15  and  16  for  N  =  180).  The  smoothing  of  the 
field  reduces  the  errors  associated  with  the  assumption  that  the 
topography  is  a  level  surface  (see  Table  4  and  Figure  17  for  a  numerical 
example).  The  considerations  of  Section  4.1  lead  to  a  feasible  hybrid 
model:  the  spherical  harmonic  model  accounts  for  remote  zone  effects; 

the  residual  topographic  model  accounts  for  high  frequency  variations  of 
the  field;  and  the  classical  integral  model  accounts  for  whatever  is  the 
residual  field  not  already  modeled  by  the  first  two  models.  The 


130 


smoothness  of  the  residual  field  allows  the  minimization  of  the  errors 
associated  with  assumption  of  the  classical  integral  model,  namely  that 
the  earth’s  surface  is  a  level  surface. 

Section  4.2  summfu'izes  the  equations  and  operational  programs 
associated  with  the  above  three  complementary  models.  Then,  numerical 
tests  are  conducted  in  Section  4.3,  over  a  mountainous  test  area  in  New 
Mexico.  Radial  and  horizontal  disturbance  values  are  computed  along 
five  vertical  test  lines,  at  altitudes  5,  10,  30,  100,  and  500  km.  The 
tests  are  designed  to  (a)  confirm  earlier  discussions  and  conclusions 
(see  Tests  1,  2,  4,  7,  8);  (b)  establish  operational  details  of 
implementation  (see  Tests  3,  5,  6);  and  (c)  assess  the  quality  of  the 
results  (see  Test  9).  An  auxiliary  test  to  Test  9  is  conducted  (Section 
4.3.4),  concluding  a  significant  improvement  caused  by  accounting  for 
the  mean  topography  in  the  classical  integral  models  (see  Table  6). 

Test  9  shows  errors  (maximal  value  =  8  mgals)  that  can  be 
encountered  if  the  classical  integral  model  is  used  to  represent  the 
residual  topographic  field.  Part  of  the  errors  is  caused  by  the 
neglection  of  the  variations  of  the  topography  with  respect  to  the  mean 
topography.  However,  as  to  be  indicated  in  the  tests  of  Chapter  5,  most 
of  the  errors  can  be  attributed  to  the  inherently  less  resolution  of  the 
data  used  in  the  integral  model  (i.e.,  5'x5'  anomalies)  as  compared  with 
the  data  used  to  generate  the  residual  topographic  model  (i.e.,  30"x30" 
topographic  heights). 

To  eliminate  the  theoretical  modeling  error  caused  by  the  neglection 
of  topographic  variations  around  the  mean  elevation,  we  can  resort  to 
the  use  of  the  so-called  discrete  approaches.  These  approaches  arise 
from  solving  the  Bjerhammar  problem.  Two  major  discrete  approaches 
are  studied  in  Chapter  5,  namely,  the  Dirac  approach  and  the  least 
squares  collocation  approach.  Three  versions  of  the  Dirac  approach  are 
experimented  on,  those  based  on  solving  on  the  internal  sphere  the  (a) 
gravity  anomaly  impulses  &g**,  (b)  point  masses  and  (c)  point  dipoles 
P*.  For  the  least  squares  collocation  (l.s.c.)  approach  two  versions  are 
experimented  on,  those  based  on  generating  the  empirical  covariance 
function  from  (a)  white  noise  gravity  anomaly  Ag*;  and  (b)  white  noise 
disturbing  potential  T*.  All  these  discrete  approaches  rigorously  take 
as  input  data  a  finite  number  of  gravity  anomalies  given  on  the  earth’s 
surface.  The  techniques  then  produce  artificial  continuity  from  the 
finite  data  by  postulating  the  form  of  the  harmonic  approximating 
function  to  T  (see  (5.17)  for  the  Dirac  case,  and  (5.65)  for  the  l.s.c. 
case)  and  fitting  this  form  to  the  given  data. 

The  Dirac  (Ag*)  and  Dirac  in*)  produced  practically  the  same  results 
in  the  numericcd  tests,  because  their  base  functions  behave  similarly  in 
the  limit  as  n  Such  high  degree  harmonics  are  the  ones  important 

for  local  modeling  problems.  The  Dirac  (^*)  required  a  deeper 
Bjerhammar  sphere  compared  with  the  Dirac  (Ag*  or  >i*)  because  of  the 
faster  decay  of  the  base  function  in  the  Dirac  (Ag*  or  /i*)  (see  (5.26), 
(5.33),  (5.37)).  Exactly  the  same  situation  exists  for  the  l.s.c.  systems. 


The  l.s.c.  (T*)  system,  which  as  noted  under  (5.92)  behaves  similarly  as 
a  l.s.c.  system,  required  twice  as  deep  Bjerhammar  sphere  as  the 

l.s.c.  (Ag*)  system  (see  (5.91),  (5.92),  Figure  26)  in  order  to  generate  the 
same  correlation  lengths  on  the  earth’s  surface. 

Comparing  the  Dirac  and  l.s.c.  systems,  we  numerically  found  the 
former  to  result  in  more  stable  linear  equation  systems  (see  also 
Sjoberg,  1978).  For  the  relatively  large  0*.5x0*.5  spacing  of  surface  data 
used  in  our  tests  in  Canada,  no  matrix  stability  problems  were 
encountered  and  all  the  Dirac  and  l.s.c.  systems  tested  produced  the 
same  quality  of  results  at  their  optimal  depths  to  the  Bjerhammar 
sphere.  For  the  dense  5'x5'  data  used  in  our  New  Mexico  tests, 
however,  the  Dirac  results  were  superior  to  those  of  the  l.s.c.  because 
of  the  ill-conditioned  matrices  in  the  latter  system. 

Considering  the  rest  of  the  tests  in  Section  5.4.3  and  summing  up 
the  report:  we  recommend  the  hybrid  model  of  Chapter  4  for 
operational  modeling  of  the  spatial  disturbance  vector,  with  possible 
replacement  of  the  classical  integral  model  by  a  collocation  procedure  for 
accurate  determinations  over  mountainous  areas.  The  residual 
topographic  model  should  always  be  used  whenever  a  detailed  DTM  is 
available,  since  both  the  integral  and  collocation  models  sire  necessarily 
limited  in  resolution  through  the  use  of  gravity  anomaly  data  with 
practical  spacing  that  can  only  be  expected  to  be  around  5'x5'.  The 
limitation  of  resolution  can  be  a  more  significant  concern  than  the 
non-rigorous  accounting  of  the  topography  in  the  integral  model  (see 
pairings  (1),  (2),  and  (9)  of  Table  11).  The  integral  model  with  mean 
topography  accounted  for  is  in  reasonably  good  agreement  with  a 
collocation  procedure  (the  Dirac  model,  see  pairing  (9)  of  Table  11). 
Matrix  conditioning  problems  with  the  l.s.c.  approach  support  preference 
to  the  Dirac  systems  for  rigorous  treatment  of  the  topography  at 
detailed  (5'x5')  resolutions. 

A  collected  summary  of  the  suggested  computational  procedures 
discussed  in  this  report,  with  references  to  pertinent  equations  found  in 
the  body  of  the  report,  is  given  as  an  Appendix. 


APPENDIX 


Summary  of  Suggested  Computational  Procedures  for  the  Computation  of 
the  Gravity  Disturbance  Components  in  Space  from  Surface  Gravity 

Anomaly  Data. 


A.  The  Spherical  Harmonic  Model  (SHM) 

1.  To  compute  spherical  harmonic  potential  coefficients  C*^  from  surface 
gravity  anomalies  dg,  apply  these  steps: 

(a)  reduce  dg  to  gravity  anomalies  dg^  on  the  ellipsoid,  using  (2.3); 

(b)  C(»pute  conventional  potential  coefficients  C„,„  froa  dg^ 
employing  ellipsoidal  corrections,  using  (2.9)  and  (2.54)-(2.56), 
with  the  RnmV  denoting  the  quantity  in  braces  in  (2.32)  and  the 
Sn„g  being  given  by  (2.51),  This  solution  implies  the  use  of  up 
to  second  radial  derivative  of  the  dg-field  to  analytically 
continue  dg^  to  gravity  anomalies  dg,  on  the  equatorial  sphere, 
and  should  be  sufficient  for  n  up  to  the  envisioned  n  -  360; 

(c)  fully  normalize  the  using  (2.19)  to  finally  obtain 

2.  See  Figures  4  and  5  to  know  approximately  how  much  information  is 

left  unresolved  by  the  available  degree  Nr.f  of  SHM  at  the  altitude 
H  of  interest.  For  example,  given  =  180  and  H  =  30  km  we 

find  a  6-mgal  RMS  signal  left  to  be  resolved  by  local  models.  Since 
the  total  RMS  signal  at  H  =  30  km  is  seen  to  be  26  mgals,  we  then 
have  6/26  =  23X  of  the  signal  energy  left  unresolved  using  = 

180. 

3.  Decide  on  how  large  cap  of  local  gravity  anomaly  data  should  be 

used  to  complement  the  SHM.  This  can  be  judged  from  the 
magnitude  of  truncation  error  that  can  be  tolerated.  The  truncation 
error  in  this  case  consists  of  the  commission  and  omission  errors 
due  to  the  use  of  the  SHM  to  represent  the  remote  zone.  The 

truncation  error  is  expressed  by  (3.52)  for  the  radial  disturbance 
and  (3.62)  for  the  total  horizontal  disturbance.  Anticipating  the  use 
of  a  high  degree  SHM  (N^.f  ^  180)  then  we  have  shown  that  it  will 
be  sufficient  to  use  unmodified  truncation  coefficients  for  the  values 
kjl  and  vj{  needed  in  (3.52)  and  (3.62).  These  unmodified  coefficients 
can  be  computed  using  the  subroutines  of  Shepperd  (1979).  A 
general  idea  of  reasonable  cap  sizes  for  given  altitudes  and 
accuracies  may  be  obtained  from  the  discussions  and  figures  of 
Chapter  3  or  from  earlier  experience,  e.g.,  in  Rapp  (1966b)  and 
Sunkel  (1981a). 
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4.  To  generate  a  dense  regular  grid  of  surface  gravity  anomalies  4g® 

implied  by  the  SHM,  use  (4.12)-(4.14>  and  the  fast  program  by  Rizos 
(1979).  The  Ag®  are  to  be  subtracted  from  the  originally  given 

surface  gravity  anomaly  data  which  should  also  have  been 
established  on  a  regular  grid  using  procedures  described  in  Cruz 
and  Laskowski  (1984,  sec.  8). 

5.  To  generate  the  spatial  disturbance  components  implied  by  the  SHM, 
use  (4.15)-(4.17).  Point  by  point  evaluations  can  be  done  using  the 
program  by  Rapp  (1982b).  Evaluations  on  a  regular  grid  in  a 
limited  area  can  be  done  efficiently,  using  a  program  that  we 
modified  from  Rizos  (1979)  to  compute  all  three  disturbance 
components  in  a  single  run. 


B.  The  Residual  Topographic  Model  (RTM) 

1.  To  generate  rigorous  values  of  surface  gravity  anomalies  Ag^  implied 
by  the  RTM,  use  the  program  of  Forsberg  (1984).  To  limit  the 
expense  in  computer  time,  an  option  exists  to  use  the  approximate 
Ag*  given  by  (4.28)  or  even  (4.29)  in  an  outer  zone,  e.g.,  outside  a 
1* -radius  from  the  spatial  points  at  which  the  disturbances  are  to 
be  computed.  Testa  directly  related  to  the  use  of  an  approximate 
Agt  are  given  as  Tests  5  and  6  in  Chapter  4.  The  Ag^,  like  the  Ag= 
of  the  SHM,  are  to  be  subtracted  from  the  surface  gravity  anomalies 
to  form  residual  gravity  anomalies  that  continue  to  refer  to  the 
earth’s  surface  (see  (4.18)). 

2.  To  generate  the  spatial  disturbance  components  implied  by  the  RTM, 
again  use  the  program  of  Forsberg  (ibid.).  Introduce  the 
transformation  (4.30)  into  Forsberg's  program  to  compute  the 
North-South  and  East-West  disturbance  components  instead  of  the 
deflections  of  the  vertical. 


C.  The  Integral  Model 

1.  The  residual  gravity  anomaly  Ag'^  of  (4.18)  can  be  input  into  the 
integral  model  to  complete  the  modeling  of  the  disturbance  vector  in 
its  entire  frequency  range.  The  Ag°  are  assumed  to  refer  to  an 
equipotential  surface,  one  passing  through  the  mean  topographic 
elevation  h„  in  the  local  area  of  computations.  The  h„  is  specifically 
used  in  (4.25). 

2.  The  actual  generation  of  the  spatial  disturbance  components  is 
performed  in  the  integral  model  using  (4.19)-(4.21).  A  program 
modified  from  Rapp  (1966b)  can  be  used  for  this  purpose.  Our 
modified  program  accepts  5'x5'  mean  anomalies,  and  the  so-called 
integrated  kernel  evaluation  was  implemented  for  those  5'x5'  blocks 
whose  centers  were  within  a  distance  of  10'  from  a  computation 
point.  The  last  procedure  was  implemented  to  avoid  large  kernel 
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discretization  errors  associated  with  the  use  of  the  center  point 
kernel  evaluation  at  low  altitudes. 

3.  The  disturbance  components  computed  from  the  integral  model  can  be 
added  to  those  from  the  spherical  harmonic  and  residual  topographic 
modelSi  to  form  the  final  disturbance  values  in  space  (see  Table  7). 

Note:  If  the  RTM  had  not  been  used  as  one  of  the  complementary 
modelsi  then  the  part  of  the  field  that  would  have  been  represented 
by  the  RTM  would  in  this  case  be  included  in  the  field  represented 
by  the  integral  model.  Errors  would  be  introduced  in  this  way,  the 
magnitudes  of  which  were  studied  over  mountainous  area  in  Test  9 
of  Chapter  4  (see  also  row  V  of  Table  7  and  pairing  (1)  of  Table  11). 
It  has  been  concluded  that  a  large  part  of  the  observed  errors  is 
due  to  the  unavoidable  loss  of  resolution  caused  by  the  use  of  5'x5' 
anomalies  in  the  integral  model  as  opposed  to  the  use  of  30"x30" 
height  information  in  the  RTM.  A  smaller  (but  still  significant)  part 
of  the  error  is  caused  by  the  non-rigorous  treatment  of  the 
topography  in  the  integral  model,  and  this  error  source  could  in 
theory  be  avoided  through  the  use  of  collocation  techniques  in 
space  (for  a  numerical  study  of  this  point  compare,  e.g.,  pairings  (1) 
and  (2)  of  Table  11). 


D.  The  Collocation  Model 

1.  Motivation  for  Use.  For  high  accuracies  in  mountainous  areas  it  is 
recommended  to  employ  collocation  techniques  to  rigorously  account 
for  the  shape  of  the  topography.  For  economy  and  improved 
convergence  of  the  solution  it  would  again  be  advisable  to  use  the 
RTM  to  remove  most  of  the  roughness  of  the  original  field,  leaving  a 
relatively  smooth  residual  field.  The  data  would  then  be  surface 
residual  anomalies  4g<=,  which  can  in  fact  be  the  of  (4.18). 

2.  Summary  of  our  Results.  We  have  seen  no  preference  of  one  over 
the  other,  among  the  use  of  gravity  anomaly  impulses,  point  masses, 
and  point  dipoles  under  the  Dirac  approach  to  collocation.  We  have 
concluded  a  preference  of  the  Dirac  over  the  least  squares 
collocation  systems,  for  reason  of  stability  of  linear  equation  system 
at  such  high  resolution  as  5'x5'.  The  base  functions  and 
covariance  functions  needed  for  disturbance  computations  under  the 
studied  collocation  systems  are  given  in  Chapter  5. 

3.  Operational  Procedures.  Either  use  collocation  to  replace  the 
integral  model  entirely,  or,  for  reasons  of  computer  economy,  one 
may  choose  to  simply  complement  the  integral  model  by  collocation  in 
an  inner  zone  close  to  the  computation  points.  For  the  latter 
option,  numerical  implementation  may  be  patterned  eifter  Lachapelle 
(1977)  but  using  the  base  function  or  covariance  function  desired. 
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Note  that  in  our  studies  we  were  specifically  concerned  with 
the  verification  and  comparison  of  various  collocation  solutions,  not 
with  the  development  of  fully  operational  procedures  for 
implementing  such  solutions.  Therefore,  we  worked  under  controlled 
conditions:  regularly  gridded  data  and  a  high  frequency  field  with 
negligible  remote  zone  effects  outside  the  small  2*x2*  test  area  (see 
Chapter  5).  In  operational  cases  we  would  like  to  efficiently  apply 
collocation  using  data  in  larger  areas  (see,  e.g.,  Blaha  (1983)), 
and/or  be  able  to  establish  and  use  variable  data  density  depending 
on  the  roughness  of  particular  sections  of  the  computation  area 
(see,  e.g.,  Sunkel  (1983a)).  Other  important  issues  are  the 
operational  establishment  of  optimal  depth  to  the  Bjerhammar  sphere, 
and  the  formation  and  numerical  solution  of  the  relevant  linear 
equation  systems.  For  developing  operational  collocation  procedures 
we  refer  the  reader  to  recent  works,  of  which  we  mention  Blaha 
(1983),  Brennecke  and  Lelgemann  (1984),  Sunkel  (1983a,  1983b), 
Tscherning  and  Forsberg  (1983). 
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